Award  Number:  W81XWH-06-1-0235 


AD 


TITLE:  Automated  Patient  Positioning  Guided  by  Cone-Beam  CT  for  Prostate 

Radiotherapy 

PRINCIPAL  INVESTIGATOR:  Ming  Chao,  Ph.D. 


CONTRACTING  ORGANIZATION:  Stanford  University 

Stanford,  CA  94305-5401 


REPORT  DATE:  January  2008 


TYPE  OF  REPORT:  Annual 


PREPARED  FOR:  U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


DISTRIBUTION  STATEMENT:  Approved  for  Public  Release; 

Distribution  Unlimited 


The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and 
should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision 
unless  so  designated  by  other  documentation. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

01-01-2008  Annual  1  JAN  2007  -  31  DEC  2007 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 


Automated  Patient  Positioning  Guided  by  Cone-Beam  CT  for  Prostate  Radiotherapy  5b.  grant  number 

W81XWH-06- 1-0235 

5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 

Ming  Chao,  Ph.D. 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


E-Mail:  mingchao@stanford.edu 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Stanford  University 
Stanford,  CA  94305-5401 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


14.  ABSTRACT 

Modem  radiotherapy  equipment  is  capable  of  delivering  high  precision  conformal  dose  distributions  to  the  target.  However,  the  target  localization  especially  for 
soft-tissue  target  such  as  prostate  is  an  issue,  because  of  its  possible  non-rigid  internal  motion  relative  to  bony  structures  or  external  landmarks.  Recently,  a 
new  technology  of  kV  cone-beam  CT  (CBCT)  has  been  integrated  onboard  with  the  linear  accelerator.  The  prostate  target  localization  could  be  potentially 
done  more  accurately  with  the  aid  of  the  on-board  CBCT  imager.  In  this  project  the  PI  group  developed  a  novel  technique  of  automatic  patient  positioning  for 
prostate  cancer  radiation  therapy.  A  few  important  milestones  have  been  achieved  toward  the  goal  of  the  project.  These  include:  (i)  Established  a  novel 
technique  to  enhance  on-board  cone-beam  CT  and  to  effectively  reduce  the  radiation  dose  incurred  in  the  scanning  process;  (ii)  Developed  a  robust 
registration  model  with  improved  metric  function;  (iii)  Developed  an  innovative  narrow  shell  technique  for  better  target  localization  and  model  the  critical 
structure;  (iv)  Established  the  method  for  auto-propagation  of  contours  from  planning  CT  to  CBCT  based  on  feature  based  spline  deformable  registration.  With 
more  evaluations,  the  automatic  patient  positioning  technique  and  relevant  research  outcomes  will  soon  be  part  of  routine  practice  at  Stanford  University 
Hospital.  The  data  and  experiences  we  accumulated  in  this  project  are  well  documented,  and  server  as  a  good  reference  for  effectively  utilizing  CBCT  imager, 
and  will  definitely  enhance  the  outlook  of  the  onboard  CBCT  in  guiding  radiation  therapy  as  a  whole. 


15.  SUBJECT  TERMS 

Prostate  Cancer 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

USAMRMC 

a.  REPORT 

U 

b.  ABSTRACT 

U 

c.  THIS  PAGE 

U 

uu 

95 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Table  of  Contents 


Introduction . 4 

Body . 4 

Key  Research  Accomplishments . 7 

Reportable  Outcomes . 8 

Conclusions . 9 

References . 9 

Appendices . 10 


3 


Introduction 


This  postdoctoral  traineeship  grant  (W81XWH-06-1-0235,  entitled  “Automated  patient 
positioning  guided  by  cone-beam  CT  (CBCT)  for  prostate  radiotherapy”)  was  awarded  to 
Dr.  Tianfang  Li,  the  original  principal  investigator  (PI)  on  January  1st,  2006.  Dr.  Tianfang 
Li  left  Stanford  University  and  accepted  a  faculty  position  at  the  University  of  Texas 
Southwestern  Medical  Center  at  Dallas  in  the  end  of  2006.  The  award  was  transferred  to 
the  current  PI,  Ming  Chao  Ph.D.  It  took  a  few  months  for  the  transition  paper  work  to 
complete.  A  no-cost-extension  was  filed  to  extend  the  period  of  the  project  for  one  year, 
namely,  through  December  31st,  2008.  This  is  the  annual  report  for  the  second  funding 
period  (January  1st,  2007  -  December  31st,  2007). 

The  goal  of  this  project  is  to  develop  a  clinically  practical  technique  for  prostate 
patient  positioning  based  on  newly  emerged  CBCT  on-board  imaging  system.  Under  the 
generous  support  from  the  U.S.  Army  Medical  Research  and  Materiel  Command 
(USAMRMC),  the  PI  has  gained  a  tremendous  amount  of  knowledge  on  prostate  cancer 
and  prostate  cancer  management.  The  support  has  also  made  it  possible  for  the  PI  to 
contribute  significantly  to  prostate  cancer  research.  A  number  of  conference  and  refereed 
publications  have  been  resulted  from  the  support.  In  this  report,  the  Pi’s  research 
activities  in  the  past  year  are  summarized  and  the  accomplishments  of  the  proposed  tasks 
in  the  Statement  of  Work  are  presented. 


Body 


I.  Introduction 

Modern  radiotherapy  equipment  is  capable  of  delivering  high  precision  conformal  dose 
distributions  to  the  target.  However,  how  to  precisely  locate  the  target  prior  to  each 
treatment  fraction,  especially  for  soft  tissue,  remains  a  problem,  due  to  possible  non-rigid 
internal  motion  relative  to  bony  structures  or  external  landmarks.  Recently,  a  new 
technology  of  kilo-voltage  cone-beam  CT  (CBCT)  has 
been  integrated  onboard  with  the  linear  accelerator 
treatment  machine.  Figure  1  shows  a  CBCT  on-board 
imager  on  a  Varian  Trilogy™  system.  Superior  to  the 
common  approach  based  on  the  two  orthogonal  images 
provided  by  the  mega-voltage  EPID,  CBCT  can  provide 
high-resolution  three-dimensional  (3D)  information  of 
the  prostate  and  other  pelvic  structures1.  Thus  the 
prostate  target  localization  could  be  potentially  done  Figure  1.  Varian  Trilogy  oncology 

more  accurately  with  the  aid  of  the  on-board  CBCT  system  with  on-board  CBCT  imager 

imager,  even  without  using  any  implanted  fiducials. 

However,  in  practice,  there  is  currently  still  a  general  lack  of  efficient  method  to  utilize 
the  3D  CBCT  images  for  accurate  patient  positioning.  The  two  major  inherent  difficulties 
in  the  task  of  prostate  localization  are:  (1)  the  inter-fraction  organ  motions  are  more 
pronounced  at  pelvis  than  other  sites,  due  to  the  involuntary  physiological  motions  of  the 
involved  organs  (such  as  the  rectum  and  bladder)2'  3  (Appendices  5  and  6);  (2)  the 
prostate  has  relatively  low  contrast  resolution  in  CT  images,  which  is  likely  to  cause 


4 


more  errors  in  target  delineation  or  localization4  (Appendices  1  and  2  ).  Extensive  efforts 
have  been  made  in  utilizing  CBCT  for  prostate  target  positioning,  for  instance,  to  position 
the  patient  using  bony  landmarks.  Since  the  prostate  gland  can  move  independently  to  the 
bony  structures  and  its  shape  may  change  from  time  to  time,  this  commonly  used 
positioning  technique  results  in  large  uncertainties  in  prostate  targeting.  Manually 
matching  the  prostate  target  is  possible  but  it  is  very  time  consuming  and  will  introduce 
significant  inter/intra-operator  uncertainties.  The  goal  of  this  project  is  to  develop 
effective  volumetric  image  guidance  techniques  to  facilitate  the  patient  positioning  for 
prostate  radiation  therapy. 

II.  Novel  strategies  of  enhancing  the  quality  of  CBCT  images 

Onboard  CBCT  imaging  plays  an  essential  role  in  the  success  of  prostate  IMRT5'7. 
Unfortunately,  the  image  quality  of  currently  available  CBCT  is  far  from  satisfactory  due  to  a 
number  of  issues  specific  to  the  CBCT  scan  and  image  reconstruction.  Factors  that  adversely 
affect  the  CBCT  performance  include  scatter  photons,  motion  artifacts,  and  excessive  radiation 
dose.  We  were  among  the  first  in  developing  method  to  utilize  the  prior  knowledge  of  the 
system  to  improve  the  resultant  CBCT  images.  Generally,  the  patient  has  already  had  a 
planning  CT  (and  possibly  one  or  more  CBCTs)  before  an  on-treatment  CBCT  is  acquired. 
These  imaging  data  can  be  utilized  as  a  priori  knowledge  to  enhance  the  image  quality  and  to 
reduce  the  radiation  dose  in  routine  CBCT  imaging  process.  This  has  been  successfully 
demonstrated  in  our  previous  studies  on  CT  and  PET  imaging8.  For  CBCT  imaging,  due  to  the 
significantly  reduced  number  of  projections  per  reconstruction,  the  quality  of  the  phase- 
resolved  CBCT  images  is  greatly  degraded  by  the  view-aliasing  artifacts.  Most  commonly, 
acquisitions  using  multiple  gantry  rotations  or  slow  gantry  rotation  can  increase  the  number  of 
projections  and  subsequently  improve  the  image  quality.  However,  the  price  to  pay  here  is  the 
significant  increase  of  scanning  time.  Our  CBCT  image  enhancement  method  effectively 
circumvent  this  problem  through  effective  utilization  of  the  information  contained  in  other 
phases  and/or  even  previously  obtained  planning  CT  (fan  beam  CT).  By  registering  the  prior 
data  to  the  CBCT  (in  either  image  space  or  projection  space),  we  were  able  to  enhance  the 
contrast-to-noise  ratio  (CNR)  of  the  CBCT  images  by  5—10  fold,  without  noticeable 
compromise  in  the  spatial  resolution  of  the  resultant  images.  The  developed  technique 
significantly  reduces  the  view-aliasing  artifacts  in  CBCT  images  and  leads  to  high  quality 
images  for  therapeutic  guidance.  This  technique  also  dramatically  decreases  the  radiation  dose 
of  CBCT  acquisition.  The  image  quality  enhancement  makes  it  possible  to  directly  use  CBCT 
for  accurate  dose  calculation  and  IMRT  planning. 

In  summary,  along  the  direction  of  Task  4  in  the  Statement  of  Work,  we  have  done 
substantial  work  on  improving  the  perfonnance  of  the  CBCT  system.  We  developed  methods 
for  enhancing  the  perfonnance  of  CBCT  and  reducing  radiation  dose  incurred  during  frequent 
CBCT  scans4' 9.  These  accomplishments  are  attached  in  the  Appendices  (items  1  and  2),  which 
conclude  the  task  4  in  the  Statement  of  Work. 

III.  Refinement  of  image  registration  technique 

In  the  previous  funding  period,  a  novel  automated  patient  positioning  technique  was 
developed  for  prostate  cancer  patients  using  defonnable  image  registration.  Compared  to 
the  conventional  approach,  the  new  technique  accounts  for  the  complex  non-rigid  motion 
of  the  soft-tissue  organs  and  balances  the  prostate  target  and  the  organs  at  risk 
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simultaneously.  The  method  was  verified  using  a  pelvis  phantom  consisting  of  bladder, 
prostate,  rectum  and  bony  structures.  The  study  was  detailed  in  the  previous  annual 
progress  report.  We  have  further  refined  the  registration  technique  and  significantly 
improved  the  computational  efficiency  of  the  technique,  thus  making  it  possible  for 
future  clinical  application.  Briefly,  our  efforts  were  focused  on  improving  the  registration 
technique  by  using  different  similarity  measures  such  as  nonnalized  cross  correlation  and 
mutual  information  (Appendices  1,  3,  and  4).  We  found  that  by  using  the  Mattes  Mutual 
Information  (MMI)  metric,  not  only  significant  amount  of  computing  time  was  saved,  but 
also  precision  of  the  registration  was  achieved10.  The  central  concept  of  mutual 
information  (MI)  is  the  calculation  of  entropy.  For  an  image  A,  the  entropy  is  defined  as 
H(A )  =  -J  pA(a) log pA(a)da  ,  where  p A(a)  (also  termed  as  marginal  probability  density 

function  (PDF))  is  the  probability  distribution  of  grey  values  (image  intensities)  which  is 
estimated  by  counting  the  number  of  times  each  grey  value  occurs  in  the  image  and 
dividing  those  numbers  by  the  total  number  of  occurrences.  Given  two  images,  A  and  B, 
their  joint  entropy  is  H(A,B)  =  -J  J  p 4B(a,b) log p AB{a,b)dadb  ,  where p 4B(a,b)  is  the 

joint  PDF  defined  by  a  ratio  between  the  number  of  grey  values  in  the  joint  histogram 
(feature  space)  of  two  images  and  the  total  entries.  The  mutual  infonnation  is  generally 
expressed  as  MI (A,  B)  =  H(A)  +  H(B )  -  H(A,  B ) .  MI  measures  the  level  of  information 
that  a  random  variable  (e.g.,  Ia(x))  can  predict  about  another  random  variable  (e.g.  Ib(x)). 
Different  from  the  conventional  MI,  where  two  separate  intensity  samples  are  drawn 
from  the  image,  the  Mattes  implementation,  MMI,  uses  only  one  set  of  intensity  to 
evaluate  both  the  marginal  and  joint  PDFs  at  discrete  positions  or  bins  that  unifonnly 
spread  within  the  dynamic  range  of  the  images11.  Entropy  values  were  computed  by 
summing  over  all  the  bins.  The  introduction  of  the  mutual  infonnation  metric  has  showed 
significant  improvement  on  the  registration  precision  and  helped  reducing  the  impact  of 
scatter  and  other  factors  inherent  in  CBCT  imaging. 

IV.  Contour  propagation  from  planning  CT  to  CBCT 

Large  interfraction  organ  motion  exists  in  radiation  therapy  of  pelvic  diseases,  such  as 
prostate  and  cervical  cancers.  Adaptive  therapy  provides  a  viable  option  to  ensure  an 
adequate  dose  coverage  of  the  tumor  target  while  sparing  the  sensitive  structures  ’  . 
Toward  the  general  goal  of  more  accurate  patient  positioning  and  clinically 
implementation  of  adaptive  therapy  of  prostate  cancer,  we  developed  an  effective 
regional  contour  mapping  technique  to  automatically  propagate  rectum  and  prostate 
contours  from  planning  CT  to  CBCT2, 3' 10' 14  (Appendices  3,  4,  5,  and  6).  This  effort  is 
based  on  the  precise  deformable  registration  between  planning  CT  and  CBCT  images. 
Different  from  disease  sites,  such  as  the  lungs  and  liver,  the  contour  mapping  here  in 
prostate  radiation  therapy  is  complicated  by  two  factors:  (i)  the  physical  one-to-one 
correspondence  may  not  exist  due  to  the  insertion  or  removal  of  some  image  contents 
within  the  rectum;  and  (ii)  reduced  signal  to  contrast  ration  of  the  CBCT  images  due  to 
increased  scatter  in  CBCT  scanning.  A  solution  customized  to  rectum  mapping  is 
proposed  to  overcome  the  above  two  limitations  and  allow  us  to  take  advantage  of  the 
regional  calculation  algorithm  yet  avoiding  the  nuisance  of  rectum/bladder  filling.  The 
implementation  of  the  contour  propagation  was  achieved  with  the  narrow  shell  technique. 
The  shell  was  constructed  outside  the  manually  delineated  contours  on  planning  CT 
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including  the  region  of  1  ~  2  cm.  The  central  idea  here  is  to  exclude  the  volume  inside  the 
rectum  because  its  content  may  change  from  day  to  day.  Prostate  contour  can  be  mapped 
similarly,  but  the  shell  spans  the  regions  both  inside  and  outside  the  manually  segmented 
prostate.  Figure  2  illustrates  the  planning  CT  and  the  template  contours  for  the  ROIs  such 
as  prostate  and  rectum.  Based  on  the  narrow  shell  technique,  various  deformable  models 
such  as  B-Spline  and  Thin  Plate  Spline  models  were  employed  to  achieve  the  goal  of 
contour  propagation15'  I6.  A  feature  based  model  using  the  Scale  Invariant  Feature 
Transformation  (SIFT)  was  also  explored  for  higher  registration  accuracy3,  17,  18 
(Appendix  5).  The  mapped  rectal  contour  is  illustrated  in  Figure  3.  This  task  was 
completed  during  the  second  funding  period  -  one  manuscript  reporting  the  use  of  narrow 
shell  based  contour  mapping  strategy  has  been  accepted  by  Physics  in  Medicine  and 
Biology  for  publication.  In  addition,  a  paper  reporting  the  use  of  SIFT  method  to  auto- 
identify  the  tissue  features  for  robust  deformable  registration  and  contour  mapping  in  the 
pelvic  region  has  been  conditionally  accepted  by  Medical  Physics  for  publication 
(Appendices  5  and  6). 


Figure  2.  Narrow  shell  representation 
of  regions  of  interest. 


Figure  3.  Mapped  rectal  contours  (in 
red).  Manual  contour  is  overlaid  (blue) 


Key  Research  Accomplishments 


•  Established  a  novel  technique  to  enhance  on-board  CBCT  and  to  effectively  reduce 
the  radiation  dose  (Task  4  of  SOW) 

•  Improved  the  registration  technique  by  investigating  various  similarity  measures 
including  both  the  nonnalized  cross  correlation  metric  and  the  Mattes  mutual 
information  metric  (Task  2  of  SOW) 

•  Developed  a  novel  technique  of  narrow  shell  representation  of  the  regions  of  interest 
better  target  localization  and  critical  structure  modeling  to  account  for  the  large  inter¬ 
fraction  (Task  3  of  SOW) 

•  Established  the  framework  of  contour  propagation  from  planning  CT  to  CBCT  based 
on  the  narrow  shell  technique  for  the  adaptive  prostate  radiation  therapy  (Task  3  of 
SOW) 
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•  M.  Chao,  L.  Xing,  “Auto-recontouring  of  CBCT  images  for  adaptive  therapy”,  Medical 
Physics  34,  2352  (2007) 

•  A.  de  la  Zerd,  M.  Chao,  B.  Armbrush,  Y.  Yang,  S.  Hancock,  C.  King,  T.  Li,  L.  Xing, 
“Adaptive  IMRT  for  improved  prostate  cancer  treatment”,  Proceedings  of  the  U.S. 
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Minds  in  Prostate  Cancer  Today  (IMPaCT)  meeting,  September  5  -8,  2007,  Atlanta, 
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•  M.  Chao,  Y.  Xie,  Q.  Le,  L.  Xing,  “Modeling  the  volumetric  change  of  head  and  neck 
tumor  in  response  to  radiation  therapy”,  International  Journal  of  Radiation  Oncology?  * 
Biology ’  *  Physics,  Volume  69,  Issue  3,  Supplement  1,  1  November  2007,  Page  S741 

•  T.  La,_M.  Chao,  L.  Xing,  Q.  Le,  “Evaluation  of  intrafraction  motion  in  head  and  neck 
cancer  during  radiotherapy”,  International  Journal  of  Radiation  Oncology >  *  Biology ?  * 
Physics,  Volume  69,  Issue  3,  Supplement  1,  1  November  2007,  Pages  S681  -  S682 

•  L.  Xing,  L.  Lee,  M.  Chao,  P.  Keall,  T.  La  and  Q.  Le,  “Clinical  implementation  of  on¬ 
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Conclusions 

In  the  past  funding  year,  the  PI  has  contributed  greatly  to  development  of  novel  prostate 
cancer  patient  positioning  technique.  He  has  also  gained  improved  knowledge  on  prostate 
cancer  and  prostate  radiation  therapy  treatment.  Various  factors  affecting  the  final 
positioning  decisions  are  addressed.  A  few  important  milestones  have  been  achieved 
toward  the  goal  of  the  project.  These  include:  (i)  Established  a  novel  technique  to 
enhance  on-board  CBCT  and  to  effectively  reduce  the  radiation  dose  incurred  in  the 
scanning  process;  (ii)  Developed  a  robust  registration  model  with  improved  metric 
function;  (iii)  Developed  an  innovative  narrow  shell  technique  for  better  target 
localization  and  model  the  critical  structure;  (iv)  Established  the  method  for  auto¬ 
propagation  of  contours  from  planning  CT  to  CBCT  based  on  feature  based  spline 
deformable  registration.  Integration  and  further  refinement  of  the  above  tools  are  in 
progress.  Together  with  the  accomplishments  as  described  in  the  first  annual  report:  (1) 
developed  and  clinically  implemented  an  automated  patient  positioning  strategy  and 
tested  with  phantom  experiments;  (2)  developed  image-to-projection  deformable 
registration  algorithm  to  improve  the  CBCT  image  quality  by  reducing  the  view-aliasing 
artifacts;  and  (3)  radiation  dose  reduction,  we  have  completed  most  tasks  as  described  in 
the  Statement  of  Work.  More  evaluations  and  testing  of  the  developed  tools  are  being 
carried  out.  The  data  and  experiences  accumulated  in  this  project  are  well  documented, 
and  serve  as  a  good  reference/guidance  for  effective  use  of  volumetric  CBCT  imaging 
technology.  The  research  improves  volumetric  image  guided  prostate  radiation  therapy 
and  enhances  the  outlook  of  the  onboard  CBCT  in  clinical  practice. 
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OPTIMIZING  4D  CONE-BEAM  CT  ACQUISITION  PROTOCOL  FOR 
EXTERNAL  BEAM  RADIOTHERAPY 


Tianfang  Li,  Ph.D.,  and  Lei  Xing,  Ph.D. 

Department  of  Radiation  Oncology,  Stanford  University  School  of  Medicine,  Stanford,  CA 

Purpose:  Four-dimensional  cone-beam  computed  tomography  (4D-CBCT)  imaging  is  sensitive  to  parameters 
such  as  gantry  rotation  speed,  number  of  gantry  rotations,  X-ray  pulse  rate,  and  tube  current,  as  well  as  a 
patient’s  breathing  pattern.  The  aim  of  this  study  is  to  optimize  the  image  acquisition  on  a  patient-specific  basis 
while  minimizing  the  scan  time  and  the  radiation  dose. 

Methods  and  Materials:  More  than  60  sets  of  4D-CBCT  images,  each  with  a  temporal  resolution  of  10  phases, 
were  acquired  using  multiple-gantry  rotation  and  slow-gantry  rotation  techniques.  The  image  quality  was 
quantified  with  a  relative  root  mean-square  error  (RE)  and  correlated  with  various  acquisition  settings; 
specifically,  varying  gantry  rotation  speed,  varying  both  the  rotation  speed  and  the  number  of  rotations,  and 
varying  both  the  rotation  speed  and  tube  current  to  keep  the  radiation  exposure  constant.  These  experiments 
were  repeated  for  three  different  respiratory  periods. 

Results:  With  similar  radiation  dose,  4D-CBCT  images  acquired  with  low  current  and  low  rotation  speed  have 
better  quality  over  images  obtained  with  high  current  and  high  rotation  speed.  In  general,  a  one-rotation 
low-speed  scan  is  superior  to  a  two-rotation  double-speed  scan,  even  though  they  provide  the  same  number  of 
projections.  Furthermore,  it  is  found  that  the  image  quality  behaves  monotonically  with  the  relative  speed  as 
defined  by  the  gantry  rotation  speed  and  the  patient  respiratory  period. 

Conclusions:  The  RE  curves  established  in  this  work  can  be  used  to  predict  the  4D-CBCT  image  quality  before 
a  scan.  This  allows  the  acquisition  protocol  to  be  optimized  individually  to  balance  the  desired  quality  with  the 
associated  scanning  time  and  patient  radiation  dose.  ©  2007  Elsevier  Inc. 

Image-guided  radiotherapy,  Four-dimensional  cone-beam  CT,  Optimization. 


INTRODUCTION 

Onboard  cone-beam  computed  tomography  (CBCT)  imag¬ 
ing  provides  a  convenient  means  for  accurate  patient  setup 
and  dose  verification  (1-8).  However,  when  used  for  thorax 
or  upper  abdomen  imaging,  motion  artifacts  appear  in  the 
reconstructed  images  because  of  intrascanning  organ  mo¬ 
tion  within  the  field  of  view  (FOV).  According  to  the 
International  Electric  Commission  recommendation,  on¬ 
board  CBCT  systems  have  a  limited  gantry  rotation  speed  of 
60  s  per  round.  A  complete  scan  therefore  consists  of 
projection  data  from  10  to  20  respiratory  cycles  of  the 
patient,  resulting  in  large  amount  of  inconsistency  in  the 
CBCT  projection  data.  Reconstruction  algorithms  based  on 
theories  of  static  object  lead  to  images  that  are  significantly 
deteriorated  by  motion  artifacts.  Because  of  the  prolonged 
scanning  time,  these  motion-induced  adverse  effects  in 


CBCT  are  much  more  severe  than  those  seen  in  conven¬ 
tional  CT  scans  (9).  The  artifacts  not  only  blur  the  image, 
but  also  inhibit  direct  use  of  CBCT  for  dose  calculation 
(10-12). 

Four-dimensional  (4D)  CBCT,  or  respiration-correlated 
CBCT,  groups  the  acquired  projection  data  into  several  bins 
as  according  to  their  respiratory  phases.  Each  phase  bin  is 
then  reconstructed  independently  to  obtain  a  volumetric 
image  corresponding  to  that  specific  phase  (13-15).  4D- 
CBCT  not  only  greatly  reduces  the  motion  artifacts  pre¬ 
sented  in  each  of  these  phase-resolved  images,  but  also 
provides  the  dynamic  information  of  the  patient  anatomy, 
which  is  absent  in  the  three-dimensional  (3D)  case,  and 
could  be  very  useful  in  future  4D  radiotherapy  (1). 

Because  the  total  amount  of  projections  acquired  during  a 
4D-CBCT  scan  is  divided  into  several  phase  groups,  the 
number  of  projections  available  for  each  image  reconstruc- 
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tion  becomes  much  less  than  in  the  regular  3D-CBCT  case. 
As  known  from  theory,  the  number  of  projections  needed  to 
avoid  undersampling  artifacts  is  inversely  proportional  to 
the  reconstructed  voxel  size  (16).  An  insufficient  number  of 
projections  will  lead  to  severe  view-aliasing  artifacts.  To 
increase  the  number  of  projections  for  each  phase,  two 
strategies  that  can  be  employed  are:  slowing  down  the 
gantry  rotation  speed  (SGR)  and  multiple  gantry  rotations 
(MGR).  In  addition  to  the  gantry  rotation  speed  and  the 
number  of  gantry  rotations,  the  quality  of  4D-CBCT  images 
is  also  sensitive  to  parameters  such  as  the  X-ray  pulse  rate, 
the  tube  current,  and  the  patient’s  breathing  pattern.  These 
parameters  constitute  a  multidimensional  space,  making 
the  optimal  4D-CBCT  protocol  intractable.  To  balance  the 
tradeoff  between  image  quality,  scan  time,  and  patient  ra¬ 
diation  dose,  a  patient-specific  4D-CBCT  acquisition  pro¬ 
tocol  is  highly  desirable.  This  issue  is  systematically  studied 
in  this  work,  with  the  goal  of  fully  understanding  the  influ¬ 
ences  of  these  various  scanning  parameters  to  optimize 
4D-CBCT  data  acquisition  protocol. 

METHODS  AND  MATERIALS 

4D-CBCT  data  acquisition  system 

An  Acuity  simulator  (Varian  Medical  Systems,  Palo  Alto,  CA) 
was  used  in  this  work  for  all  CBCT  imaging.  In  the  case  of  regular 
3D-CBCT  simulation,  a  gantry  rotation  speed  of  8°/s  is  used.  The 
X-ray  tube  operates  at  125  kVp  and  80  mA  with  a  pulse  width  at 
each  projection  angle  of  9  ms  when  no  bow-tie  filter  is  used.  The 
pulse  width  is  increased  to  25  ms  if  a  bow-tie  filter  is  used  to 
maintain  similar  photon  statistics.  The  data  are  acquired  at  ~15 
frames/s  and  a  full  rotation  (slightly  more  than  360°)  consists  of 
about  685  projections  corresponding  to  an  angle  interval  of  about 
0.5°.  The  radiation  dose  of  a  one-rotation  CBCT  scan  at  isocenter 
is  approximately  4.0  cGy.  The  dimension  of  each  acquired  pro¬ 
jection  image  is  397.3  mm  X  298.0  mm,  containing  1,024  X  768 
pixels.  With  a  source-to-axis  distance  of  100  cm  and  source-to- 
detector  distance  of  150  cm,  the  field  of  view  is  around  25  cm  in 
diameter  in  the  transverse  plane  for  a  full-fan  mode.  This  can  be 
close  to  50  cm  if  a  half-fan  mode  is  employed  (by  shifting  the 
flat-panel  detector  laterally).  Here  it  should  be  noted  that  the 
half-fan  mode  generally  results  in  inferior  image  quality  when 
compared  with  that  obtained  with  full-fan  mode  if  all  other  acqui¬ 
sition  parameters  are  maintained  the  same,  and  is  a  consequence  of 
the  loss  of  the  data  redundancy  as  is  in  the  full-fan  mode.  How¬ 
ever,  for  regions  such  as  the  thorax,  a  full-fan  mode  may  result  in 
truncated  images,  which  are  not  appropriate  for  dose  calculation. 
Thus  the  half-fan  mode  is  preferred  when  a  large  field  of  view  is 
required,  and  this  mode  is  therefore  used  in  this  work. 

To  generate  the  4D  image  sets,  the  acquired  CBCT  projections 
must  be  sorted  into  a  number  of  bins  according  to  their  specific 
respiratory  phases.  The  phase  information  of  each  projection  can 
be  obtained  with  the  aid  of  a  motion  tracking  system,  for  example, 
the  Real-time  Position  Management  system  (Varian  Medical  Sys¬ 
tems)  as  done  in  4D-CT  acquisition  with  conventional  diagnostic 
CT  scanners  (17-25),  or  alternatively,  it  can  be  achieved  by 
analyzing  the  acquired  raw  data  in  the  projection  space  (26).  In  the 
latter,  one  or  more  CT-opaque  fiducials  are  placed  on  the  patient 
skin  and  the  locations  of  the  fiducials  in  the  projections  are 
detected  by  a  computer  searching  algorithm.  The  coordinate  of  the 
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Fig.  1.  An  example  of  the  projection-space  phase  tracking  used  in 
this  work.  The  black  solid  squares  and  the  red  open  circles  repre¬ 
sent  the  positions  of  the  fiducial  marker  in  the  real  world  coordi¬ 
nate  system  and  projection  space,  respectively.  The  sinusoidal 
curve  obtained  from  cone-beam  computed  tomography  projections 
can  be  used  to  determine  the  phase  of  each  projection. 


fiducial  is  recorded  as  a  function  of  the  gantry  angle,  revealing  the 
motion  status  of  the  object  for  each  projection.  In  this  work,  the 
projection-space  approach  is  used.  Figure  1  shows  an  example  of 
such  a  plot  for  the  periodically  moving  phantom.  The  projections 
are  phase-tagged  according  to  the  resultant  sinusoidal  curve,  and 
sorted  into  10  phase  groups.  The  data  in  each  phase  bin  are  then 
reconstructed  independently  using  a  Feldkamp  algorithm  with  a 
pixel  size  of  0.5  mm  X  0.5  mm  in  the  cross-section  plane  and  a 
slice  thickness  of  1.0  mm. 

Motion  phantom 

To  investigate  the  influence  of  various  scanning  parameters  on 
the  quality  of  4D-CBCT,  a  motion  phantom  was  constructed.  The 
motion  phantom  consisted  of  a  commercial  CT  calibration  phan¬ 
tom  CatPhan  600  (The  Phantom  Laboratory,  Inc.,  Salem,  NY)  that 
is  placed  on  top  of  a  platform  capable  of  sinusoidal  motion  along 
three  directions:  the  phantom  moved  with  the  maximal  displace¬ 
ments  5.5  cm  in  superoinferior  direction,  1.5  cm  in  anteroposterior 
direction,  and  0.2  cm  in  lateral  direction.  The  period  of  the  motion 
was  continuously  adjustable  in  the  range  of  0.5  s-1  min.  In  this 
work,  three  different  periods  were  used  to  study  the  effects  of 
breathing  cycle  on  the  image  quality,  which  were  3.37  s,  4.59  s, 
and  6.53  s.  4D-CBCT  images  of  the  motion  phantom  were  ac¬ 
quired  using  SGR  or  MGR  methods  as  described  in  the  following 
section.  All  4D  data  acquisitions  were  using  X-ray  tube  current  of 
10  mA  unless  otherwise  stated.  For  comparison,  three  additional 
scans,  with  the  phantom  “frozen”  at  the  “peak-inspiration,”  “mid¬ 
expiration”  and  “peak-expiration”  phases,  were  obtained  as  well 
using  the  standard  3D-CBCT  protocol  (gantry  rotation  speed  at 
8°/s,  X-ray  tube  current  and  voltage  at  80  mA  and  125  kVp, 
respectively).  The  three  scans  served  as  our  “standards”  for  these 
phases. 

4D-CBCT  image  analysis 

Because  the  images  in  this  work  involve  significant  artifacts,  a 
common  metric  such  as  contrast-to-noise  ratio  (CNR),  defined  as 
CNR  =  I  S[  .S3  N  l/cr,  may  not  be  suitable,  in  the  sense  that 
a  high  CNR  output  may  not  represent  a  better  image.  For  example, 
dark  streak  artifacts  may  accidentally  decrease  the  mean  value  of 
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S2,  leading  to  a  high  CNR  value.  Thus  the  CNR  metric  may  highly 
depend  on  the  selected  regions  of  interest.  To  have  a  more  robust 
and  accurate  assessment  of  the  images  obtained  using  different 
CBCT  settings,  a  “relative  root  mean  square  error”  (denoted  by 
RE)  is  chosen  as  the  figure  of  merit  of  the  image  quality,  which  is 
defined  as 

\  1/2 

S&i)j  ,  (1) 

where  S4D  denotes  the  4D  single-phase  image,  and  S0  is  the 
standard  80-mA  3D-CBCT  images  of  the  phantom  frozen  at  the 
same  phase.  The  summation  runs  over  all  voxels  of  the  images.  In 
Eq.  1,  S0  image  is  used  as  the  “gold  standard,”  and  the  mean  square 
error  between  the  4D  images  and  the  gold  standard  is  normalized 
to  the  mean  square  of  the  true  intensity. 

Slow-gantry  rotation  strategy 

A  sufficient  number  of  projections  must  be  collected  for  each 
phase  to  remove  the  view-aliasing  artifacts  in  4D-CBCT.  One  way 
to  achieve  this  is  to  slow  down  the  gantry  rotation.  In  this  way,  the 
number  of  breathing  cycles  covered  by  a  full  rotation  is  increased, 
leading  to  more  projections  in  each  phase  group.  Figure  2  is  a 
sketch  illustrating  the  relationship  between  the  gantry  rotation 
schemes  and  the  breathing.  The  dark  areas  represent  the  available 
projections  that  belong  to  the  same  phase  group. 

To  examine  the  influence  of  gantry  rotation  speed  on  the  result¬ 
ant  image  quality,  scans  were  made  of  the  motion  phantom  for 
eight  gantry  rotation  speeds  ranging  from  8°/s  to  l°/s.  All  other 
system  variables  were  kept  constant  with  the  X-ray  tube  current 
and  voltage  in  these  eight  scans  being  held  at  10  niA  and  125  kVp, 
respectively.  This  increased  the  total  number  of  projections  from 
about  680  for  8°/s  speed  to  about  5,440  for  l°/s  speed.  Conse¬ 
quently,  the  number  of  projections  in  each  phase  bin  increased  by 
~8  times.  The  resultant  4D-CBCT  images  were  then  analyzed 
using  the  RE  metric  as  defined  in  Eq.  1. 

Another  important  quantity  for  SGR  4D  acquisition  is  the  rela¬ 
tion  between  image  quality  and  the  associated  radiation  exposure. 
For  a  constant  X-ray  pulse  width,  the  radiation  dose  is  directly 
proportional  to  the  number  of  projections  (beam-on  time)  and 
essentially  linear  to  the  X-ray  tube  current.  Hence  a  scan  with  a 
rotational  speed  of  l°/s  at  10  mA  current  will  result  in  a  similar 
amount  of  radiation  as  a  scan  with  a  2°/s  rotation  speed  at  doubled 
the  tube  current.  To  investigate  this  image  quality-tube  current 
relationship,  4D-CBCT  images  of  the  motion  phantom  were  ac¬ 
quired  for  l°/s-10  mA,  2°/s-20  mA,  4°/s-40  mA,  57s-50  mA,  and 
8°/s-80  mA,  for  each  of  the  three  motion  periods  mentioned 
previously.  These  scanning  schemes  deliver  a  similar  radiation 
dose  to  the  phantom.  A  quantitative  comparison  of  the  4D  images 
was  carried  out. 

Multiple-gantry  rotation  strategy 

An  alternative  to  SGR  is  MGR  4D  acquisition,  in  which  the 
gantry  now  rotates  multiple  times  back  and  forth.  As  with  SGR, 
the  projections  of  the  same  phase  are  again  collected  and  used  for 
4D-CBCT  image  reconstruction.  The  MGR  acquisition  scheme  in 
relation  to  the  breathing  pattern  is  sketched  in  Fig.  2c.  As  shown, 
both  SGR  and  MGR  strategies  lead  to  an  increase  in  the  number  of 
projections  over  conventional  CBCT  (Figs.  2b,  2c);  however,  the 
angular  distribution  of  the  projections  of  a  given  phase  group  is 
very  different.  In  SGR,  the  projections  belonging  to  a  particular 


breathing  phase  are  regularly  distributed  in  the  angular  space, 
whereas  in  MGR,  two  or  more  rotations  lead  to  the  projection 
becoming  irregular.  As  a  result,  now  there  is  a  chance  that,  during 
the  multiple  gantry  rotations,  projections  at  the  same  phase  and 
gantry  angle  may  occur  more  than  once,  leading  to  redundancy  of 
the  data  and  resulting  in  less  effective  4D-CBCT  reconstruction. 
The  angular  overlap  or  partial  overlap  of  the  projection  data  can, 
in  principle,  be  avoided/reduced  by  properly  selecting  the  initial 
phase  of  the  two  (or  multiple)  scans.  Generally,  as  illustrated  in 
Fig.  2d,  the  chance  for  two  successive  scans  to  overlap  is  minimal 
when  the  gantry  moves  in  opposite  directions.  One  can  imagine 
that  if  two  scans  are  in  the  same  direction,  at  every  angle,  the  two 
projections  from  the  two  scans  may  be  at  same  phase  in  the 
extreme  case.  This  scenario  will  never  happen  if  the  two  scans  are 
in  opposite  directions.  Two-rotation  MGR  scans  with  the  gantry 
“back  and  forth”  were  performed  for  each  of  the  previously  men¬ 
tioned  three  motion  periods  and  the  results  were  compared  with 
that  obtained  using  SGR.  To  maintain  a  constant  radiation  expo¬ 
sure,  the  gantry  rotation  speed  in  the  MGR  approach  was  increased 
by  a  factor  of  two  as  compared  with  that  of  a  single  rotation  scan. 

Relative  gantry  speed 

Four-dimensional-CBCT  image  quality  depends  not  only  on  the 
gantry  rotation  speed,  but  also  on  the  patient  respiratory  pattern. 
This  suggests  that  a  patient-specific  setting  might  be  needed  to 
achieve  a  similar  quality  of  4D-CBCT  images.  Here  we  propose  a 
new  measure  of  relative  speed  to  combine  the  two  variables  of 
gantry  rotation  and  patient  respiration,  which  is  defined  as  follows 

Relative  Speed  =  Respiratory  Period(s) 

X  Gantry  Rotation  Speed(deg/s).  (2) 

Relative  speed  has  a  unit  of  degree,  which  is  the  spanned  angle 
of  the  CBCT  scan  during  the  time  of  one  respiratory  cycle. 

The  data  acquired  in  these  studies,  for  different  motion  periods 
and  different  scanning  speed  at  constant  X-ray  tube  current  (10 
mA)  and  voltage  (125  kVp),  were  reanalyzed  to  illustrate  the 
utility  of  the  new  concept  in  optimizing  the  4D-CBCT  image 
acquisition  protocol  on  a  patient-specific  basis. 

RESULTS 

4D-CBCT  image  quality  as  a  function  of  gantry 
rotation  speed 

Four-dimensional-CBCT  images  acquired  with  eight  dif¬ 
ferent  gantry  rotation  speeds  and  constant  X-ray  tube  cur¬ 
rent  (10  mA)  for  the  phantom  moving  at  a  period  of  3.37  s 
are  summarized  in  Fig.  3.  From  the  top-left  to  the  bottom- 
middle,  the  images  correspond  to  gantry  rotation  speeds 
from  8°/s  to  l°/s.  Without  loss  of  generality,  only  the  phase 
0%  images  of  the  4D-CBCT  acquisitions  are  compared 
here.  Note  that  all  other  imaging  parameters  such  as  the  tube 
current  and  voltage  were  the  same  for  all  the  eight  scans. 
The  golden  standard  image  (obtained  with  a  standard  3D- 
CBCT  protocol  and  without  motion)  is  also  shown  at  the 
bottom  right  in  Fig.  3.  As  expected,  it  is  found  that  the 
image  quality  improves  with  decreasing  gantry  rotation 
speeds.  Similar  trends  were  also  observed  for  the  phantom 
moving  with  the  other  two  periods  (4.93  s  and  6.53  s). 
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Fig.  2.  Illustration  of  four-dimensional  cone-beam  computed  tomography  acquisitions  with  gantry  rotation  speed  (SGR) 
and  multiple  gantry  rotations  (MGR)  methods.  The  shadow  in  the  diagrams  represents  the  projections  available  for 
reconstruction  of  one  particular  phase,  (a)  Single  rotation  with  rotation  speed  of  a  conventional  three-dimensional  scan, 
in  which  the  number  of  projections  in  a  phase  bin  is  relatively  small;  (b)  with  SGR,  a  scan  expands  more  respiratory 
cycles,  hence  more  projections  of  the  same  phase  can  be  collected;  (c)  with  MGR,  the  number  of  projections  in  a  phase 
group  increases  as  shown  by  the  gray  and  blue  shadow;  (d)  the  phase  distribution  in  two  successive  rotations  of  the  MGR 
method.  Note  that  if  two  projections  overlap  at  a  particular  angle,  for  example,  for  0%  phase,  the  projections  at  neighbor 
angles  will  be  in  different  phases  because  the  two  rotations  are  in  opposite  directions. 


The  quantitative  relation  between  image  quality  and  gan¬ 
try  rotation  speed  was  analyzed  using  the  image  metric  of 
REs  and  the  results  were  plotted  in  Fig.  4.  For  all  three 


respiratory  periods,  the  relative  errors  are  found  to  increase 
with  increasing  gantry  rotation  speeds,  indicating  a  drop  in 
the  image  quality.  Also  from  Fig.  4,  it  is  seen  that  at  a  given 
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Fig.  3.  From  the  top-left  to  the  bottom-right  are  the  four-dimen¬ 
sional  cone-beam  computed  tomography  (4D-CBCT)  images  of 
phase  0%  obtained  by  varying  the  gantry  rotation  speed  from  8°/s 
to  l°/s,  and  the  three-dimensional  CBCT  images  of  the  phantom 
“frozen”  at  the  same  phase.  All  4D-CBCT  were  acquired  with 
X-ray  tube  current  of  10  mA. 


gantry  rotation  speed,  the  image  quality  improves  (smaller 
REs)  with  increasing  phantom  velocities.  This  is  not  sur¬ 
prising  because  a  full  scan  at  fast  phantom  speed  will 
contain  more  respiratory  cycles,  leading  to  a  reduction  of 
the  gap  (Fig.  2)  between  the  projections  of  same  phase  in 
two  successive  cycles,  hence  less  view-aliasing  artifacts.  To 
help  visualize  this  effect.  Fig.  5  shows  the  0%  phase  images 
obtained  for  the  three  different  motion  modes  with  the  same 
gantry  rotation  speed  of  2°/s. 
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Fig.  4.  Relation  between  relative  root  mean-square  error  and 
gantry  rotation  speed  for  three  motion  modes  with  periods  of 
3.37  s,  4.93  s,  and  6.53  s,  respectively. 


Fig.  5.  Four-dimensional  cone-beam  computed  tomography  im¬ 
ages  of  phase  0%  obtained  at  the  same  gantry  rotation  speed  of 
2°/s,  but  different  phantom  motion  periods. 


SGR  acquisition  under  the  condition  of  a  constant 
radiation  exposure 

As  mentioned  previously,  the  photon  flux,  noise  level  and 
patient  radiation  exposure  are  all  related  to  the  tube  current. 
This  is  investigated  by  comparison  of  the  0%  phase  4D 
images  for  scans  of  8°/s  at  80  mA,  5°/s  at  50  mA,  4°/s  at  40 
mA,  2°/s  at  20  mA,  and  l°/s  at  10  mA.  Because  the 
radiation  exposure  is  approximately  proportional  to  the  tube 
current,  scans  with  these  eight  settings  deliver  a  similar 
level  of  radiation  dose  to  the  phantom.  Again,  the  image 
metric  of  RE  was  employed  to  quantify  the  difference 
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Gantry  Rotation  Speed  (deg/sec) 

Fig.  6.  The  effects  of  increasing  tube  current  and  gantry  rotation 
speed.  With  the  same  radiation  dose,  lower  speed  results  in  higher 
quality  four-dimensional  images  as  a  result  of  increased  number  of 
projections  per  phase. 


among  the  images.  As  shown  in  Fig.  6,  the  image  acquired 
with  the  combination  of  l°/s  and  10  mA  setting  had  the  best 
quality  (smallest  RE).  For  visual  comparison,  the  l°/s — 10 
mA  and  8°/s — 80  mA  images  are  shown  in  Fig.  7.  Moving 
to  slower  gantry  rotation  speed  will  lead  to  each  phase 
accumulating  more  projection  data  that  can  be  used  for 
reconstruction  as  compared  with  faster  scans.  As  a  conse¬ 
quence,  the  integral  photon  number  (filtered  back  projec¬ 
tion)  at  a  spatial  point  is  not  reduced  in  the  SGR  low-current 
acquisition  as  compared  with  high-current,  high-speed 
scans.  Furthermore,  at  low  scan  speed,  the  data  become 
more  evenly  spared  along  the  360°  circle  for  each  phase, 
leading  to  less  view-aliasing  artifact  and  consequently  to 
higher  image  quality. 

As  shown  in  Fig.  6,  it  is  found  that  RE  increases  almost 
linearly  as  a  function  of  the  scan  speed.  Because  the  total 
number  of  projections  in  a  scan,  and  therefore  the  number  of 
projections  in  each  phase,  varies  linearly  with  the  gantry 
rotation  speed,  it  seems  that  the  RE  measure  is  an  appro¬ 
priate  image  metric,  which  has  a  linear  relationship  to  the 
number  of  projections  available  for  generating  the  phase 
image. 

4D-CBCT  imaging  with  single-  and 
double-gantry  rotations 

Both  SGR  and  MGR  resulted  in  an  increased  number  of 
projections.  However,  as  illustrated  in  Figs.  2b  and  2c,  the 
angular  distribution  of  the  projections  is  different  for  the  two 
acquisition  schemes.  In  general,  the  efficiency  of  the  MGR 
method  depends  on  the  level  of  projection  redundancy  or 
overlap  of  the  data  from  the  two  or  more  gantry  rotations. 
The  4D-CBCT  images  acquired  with  SGR  and  MGR  were 
compared  using  the  RE  metric,  and  the  results  are  summa¬ 
rized  in  Fig.  8.  Here  it  is  found  that  for  all  three  motion 
modes  (3.37  s,  4.93  s,  and  6.53  s,  respectively)  and  for 
different  combinations  of  the  gantry  rotations  speed  and 
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number  of  rotations,  SGR  almost  always  led  to  a  lower  RE 
or  better  image  quality  than  the  MGR  method. 

4D-CBCT  image  quality  as  a  function  of  the 
relative  speed 

The  relationship  between  the  image  quality  and  the  rela¬ 
tive  speed  is  presented  in  Fig.  9.  The  data  shown  in  Fig.  9 
were  fitted  with  a  second  order  polynomial.  The  parameters 
of  the  model  and  the  goodness  of  the  fit  are  also  shown  in 
the  figure.  The  mono  tonic  dependence  of  the  RE  on  the 
relative  speed  suggests  that  the  relative  speed  is  a  meaning¬ 
ful  concept  in  characterizing  the  4D-CBCT  data  acquisi¬ 
tion.  The  quantity  “condenses”  two  important  parameters 
(patient-specific  respiration  period  and  the  gantry  rotation 
speed)  into  one  and  is  useful  for  us  to  optimize  the  4D- 
CBCT  acquisition  protocol.  In  reality,  because  the  respira¬ 
tory  period  of  a  patient  can  be  measured  with  the  aid  of  a 
real-time  position  management  signal  or  similar  device 
( e.g .,  strain  gauge  signal)  before  the  4D-CBCT  scanning,  a 
suitable  gantry  speed  can  be  derived  to  achieve  a  prespeci¬ 
fied  image  quality  before  4D-CBCT  scan.  Alternatively,  one 
can  use  the  curve  as  shown  in  Fig.  9  to  preestimate  the  4D 
image  quality  for  an  acquisition  setting  and  to  estimate  the 
corresponding  patient  radiation  dose. 

DISCUSSION 

Cone-beam  CT  volumetric  imaging  integrated  with  a 
medical  linear  accelerator  opens  new  avenues  for  improving 
current  radiation  oncology  practice.  In  reality,  CBCT  has 
two  important  applications:  patient  setup  and  dose  recon¬ 
struction/verification.  By  imaging  the  patient  routinely  dur¬ 
ing  a  course  of  radiation  therapy,  the  accuracy  of  the  patient 
setup  can  potentially  be  improved.  Furthermore,  the  CBCT 
provides  a  pretreatment  patient  model  on  which  the  dose 
calculation  can  be  performed  using  intended  fluence  maps 
from  the  planning  system  or  other  means.  Both  applications 
rely  on  the  fidelity  and  quality  of  the  volumetric  images. 
When  imaging  in  the  region  of  thorax  and  upper  abdomen, 
however,  respiration-induced  artifacts,  such  as  blurring, 
doubling,  streaking,  and  distortion  are  observed,  which 
heavily  degrade  the  image  quality,  and  affect  the  target 


Fig.  7.  Comparison  of  four-dimensional  cone-beam  computed  to¬ 
mography  images  acquired  with  speed  l°/s  tube  current  10  mA 
(left)  and  speed  8°/s  tube  current  80  mA  (right). 
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Fig.  8.  The  effects  of  increasing  the  number  of  rotations  and  gantry  speed.  With  the  same  radiation  dose,  for  most  of 
the  cases,  low  rotation  speed  results  in  a  better  image  quality  than  multiple  rotations. 


localization  ability  and  the  accuracy  of  dose  verification. 
These  artifacts  in  CBCT  are  much  more  severe  than  those 
found  in  conventional  CT  exams  because  of  the  relatively 
slow  gantry  rotation.  In  conventional  CT,  each  rotation  of 
the  scan  can  be  completed  within  0.5  s,  during  this  period 
the  organ/tumor  motion  is  relatively  small.  Patient  body 
restraints  and  breath-hold  techniques  can  be  used  to  mini¬ 
mize  the  motion  if  necessary.  On  the  contrary,  in  CBCT 
scan,  the  gantry  rotation  speed  is  much  slower,  typically 
around  1  min  for  a  full  360°  scan  in  acquiring  the  projection 
data,  which  covers  more  than  10  breathing  cycles  for  most 
patients.  4D-CBCT  or  phase-correlated  CBCT  is  an  effec¬ 
tive  way  to  reduce/eliminate  the  motion  artifacts  and  makes 
it  useful  for  guiding  the  patient  setup  and  dose  validation  in 
the  presence  of  organ  motion. 

Radiation  dose  is  an  important  issue  in  CT  imaging  and, 
in  particular  the  onboard  CBCT  imaging  because  of  the 
repeated  use  of  modality  for  a  given  patient.  We  found  that, 
as  we  increase  the  number  of  projections  per  phase  by 
slowing  down  the  gantry  rotation  speed  or  multiple  gantry 
rotations,  the  tube  current  can  be  lowered  accordingly.  In 
this  way,  the  4D  image  quality  is  generally  not  compro¬ 


mised  and  one  can  obtain  decent  4D-CT  images  without 
increasing  the  radiation  exposure.  In  a  sense,  this  scheme  is 
to  “spread"  the  photons  of  a  projection  in  3D-CBCT  to  a 
range  of  angles  in  4D-CBCT  imaging,  which  represents  a 
better  tradeoff  between  the  signal-to-noise  ratio  and  the 
reduction  in  motion  artifacts.  In  addition  to  the  reduced 
patient  dose,  low  current  is  essential  to  avoid  overheating 
the  X-ray  tube,  especially  for  an  onboard  imager  without  oil 
cooling  system.  Although  the  4D  image  quality  directly 
relates  to  the  number  of  projections  available  for  each 
phase,  there  is  a  limit  for  slowing  down  gantry  speed. 
Beyond  this  limit,  the  time  needed  to  complete  the  scan  may 
be  too  long  to  cause  patient  discomfort  and  extra  intrascan 
motion  (other  than  the  respiration-induced  motion)  may 
lead  to  undesirable  artifacts.  There  is  a  practical  need  to 
keep  the  scanning  time  as  short  as  possible  while  maintain¬ 
ing  an  acceptable  image  quality.  The  result  shown  in  Fig.  9 
provides  guidance  for  the  4D-CBCT  acquisition  on  a  pa¬ 
tient-specific  basis,  which  predicts  the  quality  of  the  4D 
images  for  a  certain  gantry  rotation  speed  and  patient  respi¬ 
ratory  period.  Though  patient  irregular  breathing  patterns 
may  alter  the  number  and  angular  distribution  of  the  pro- 


1218 


I.  J.  Radiation  Oncology  •  Biology  •  Physics 


Fig.  9.  Relationship  between  image  quality  and  relative  speed  of 
the  four-dimensional  cone-beam  computed  tomography  scan.  The 
data  were  acquired  for  three  motion  modes  with  gantry  rotation 
speed  varying  from  l°/s  to  87s  and  constant  X-ray  tube  current 
(10  mA)  and  voltage  (125  kVp). 

jections  in  each  phase  and  thus  change  the  image  qualities, 
Fig.  9  will  nonetheless  represent  an  approximate  relation 
between  the  image  quality  and  the  average  patient  respira¬ 
tory  period.  Obviously,  this  approximation  is  valid  only 
when  the  irregularity  of  the  breathing  pattern  is  not  far  from 
ideal  periodic  motion.  It  is  also  worth  noting  that  the  influ¬ 
ence  of  breathing  irregularity  on  the  4D  images  is  different 
from  that  in  conventional  fan  beam  CT,  in  which  the  bin¬ 
ning  artifacts  may  arise  due  to  mismatch  of  the  slices 
corresponding  to  different  couch  positions. 

Although  4D-CBCT  technique  (SGR  or  MGR)  with 
phase  binning  before  reconstruction  reduces  the  motion 
artifacts  and  improves  the  image  quality,  it  deals  with  each 
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phase  independently  and  ignores  the  correlation  of  the  pro¬ 
jections  among  different  phase  groups.  It  is  of  great  interest 
to  reconstruct  a  4D  object  (with  spatial  and  temporal  pa¬ 
rameters)  rather  than  a  series  of  independent  3D  objects 
(with  only  spatial  parameters),  by  using  simultaneously  all 
the  projection  data  collected  in  different  respiratory  phase 
groups.  With  the  aid  of  deformable  registration  and  novel 
reconstruction  techniques  (9,  22,  27-33),  it  is  possible  to 
incorporate  all  projection  data  of  different  phases  into  one 
image  reconstruction  process.  In  this  way,  information  from 
the  measurement  (CBCT  scan)  is  fully  used  in  the  data 
processing,  and  more  accurate  imaging  and  better  dose 
efficiency  may  be  achieved. 

CONCLUSIONS 

We  have  demonstrated  that  4D-CBCT  images  with  sig¬ 
nificantly  reduced  artifacts  can  be  obtained  with  a  commer¬ 
cially  available  on-board  imager  system.  A  measure  of 
relative  error  was  used  to  quantify  the  image  quality.  At  the 
same  radiation  dose,  the  low-current,  low-speed  acquisition 
is  superior  to  acquisition  with  a  high  current  and  high  speed 
(<8°/s),  and,  in  most  cases,  for  the  same  total  scanning 
time,  one-rotation  low-speed  scan  is  more  preferable  as 
compared  with  an  acquisition  based  on  multiple  rotations 
and  higher  speed.  Furthermore,  we  suggested  that  4D- 
CBCT  gantry  rotation  speed  should  be  determined  on  a 
patient-specific  basis  and  have  established  a  quantitative 
relation  between  the  4D-CBCT  image  quality  and  the  rela¬ 
tive  speed.  The  approach  optimizes  the  4D  acquisition  and 
makes  an  efficient  use  of  the  available  4D-CBCT  technol¬ 
ogy.  The  study  lays  foundation  for  clinical  application  of 
4D-CBCT  and  should  have  significant  implication  to  image- 
guided  radiation  therapy. 
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Four-dimensional  (4D)  cone-beam  CT  (CBCT)  is  commonly  obtained  by  respiratory  phase  binning 
of  the  projections,  followed  by  independent  reconstructions  of  the  rebinned  data  in  each  phase  bin. 
Due  to  the  significantly  reduced  number  of  projections  per  reconstruction,  the  quality  of  the 
4DCBCT  images  is  often  degraded  by  view-aliasing  artifacts  easily  seen  in  the  axial  view.  Acqui¬ 
sitions  using  multiple  gantry  rotations  or  slow  gantry  rotation  can  increase  the  number  of  projec¬ 
tions  and  substantially  improve  the  4D  images.  However,  the  extra  cost  of  the  scan  time  may  set 
fundamental  limits  to  their  applications  in  clinics.  Improving  the  trade-off  between  image  quality 
and  scan  time  is  the  key  to  making  4D  onboard  imaging  practical  and  more  useful.  In  this  article, 
we  present  a  novel  technique  toward  high-quality  4DCBCT  imaging  without  prolonging  the  acqui¬ 
sition  time,  referred  to  as  the  “enhanced  4DCBCT”.  The  method  correlates  the  data  in  different 
phase  bins  and  integrates  the  internal  motion  into  the  4DCBCT  image  formulation.  Several  strate¬ 
gies  of  the  motion  derivation  are  discussed,  and  the  resultant  images  are  assessed  with  numerical 
simulations  as  well  as  a  clinical  case.  ©  2007  American  Association  of  Physicists  in  Medicine. 
[DOI:  10.1118/1.2767144] 

Key  words:  cone-beam,  4D  CT,  on-board  imager,  IGRT,  organ  motion 


I.  INTRODUCTION 


Medical  linear  accelerators  equipped  with  cone-beam  CT 
(CBCT)  imaging  system,  using  either  a  kV  or  MV  x-ray 
source,  have  recently  become  available  in  clinics.  The  volu¬ 
metric  image  provided  by  CBCT  opens  new  avenues  for  pa¬ 
tient  setup,  dose  verification,  and  online  treatment 
planning.  These  applications  rely  highly  upon  the  fidelity 
and  quality  of  the  CBCT  images.  When  imaging  in  the  re¬ 
gion  of  thorax  and  upper  abdomen,  however,  respiration  in¬ 
duced  artifacts,  such  as  blurring,  doubling,  streaking,  and 
distortion  are  observed,  which  heavily  degrade  the  image 
quality  and  affect  the  target  localization  ability  and  the  accu¬ 
racy  of  the  dose  calculation.1'1' 

A  recently  proposed  method  to  account  for  respiratory 
motion  during  CBCT  imaging  is  called  “respiration  corre¬ 
lated  CBCT”  or  four-dimensional  (4D)  CBCT.15-18  In  this 
method,  the  CBCT  projections  are  divided  into  several 
groups  according  to  their  respiratory  phases.  Each  phase 
group  is  then  reconstructed  independently  to  yield  a  volu¬ 
metric  image  corresponding  to  that  specific  phase.  Since  the 
selected  projections  in  each  group  have  almost  the  same 
phase,  the  method  greatly  reduces  the  motion-induced  incon¬ 
sistency  among  the  data,  leading  to  improved  image  recon¬ 
struction.  However,  the  number  of  projections  available  for 
each  image  reconstruction  is  substantially  decreased  com¬ 
pared  to  a  conventional  3D  CBCT  due  to  the  data  dividing. 
This  may  lead  to  a  serious  undersampling  problem  and  result 
in  strong  view-aliasing  artifacts  in  the  phase-resolved  im¬ 


ages.  The  artifacts  are  generated  during  the  backprojection 
step  in  the  CBCT  image  reconstruction  process  and  are  seen 
mainly  in  the  2D  axial  planes. 

To  eliminate  the  view-aliasing  artifacts  in  4DCBCT,  we 
have  recently  investigated  acquisition  techniques  of  “mul¬ 
tiple  gantry  rotation”  and  “slow  gantry  rotation”  with  low 
x-ray  tube  current.16'17  These  methods  increase  the  number 
of  projections  of  each  phase  and  significantly  improve  the 
4DCBCT  image  quality.  However  a  practical  issue  is  the 
prolonged  acquisition  time,  which  may  limit  their  clinical 
applications. 

The  purpose  of  this  work  is  to  investigate  a  novel  ap¬ 
proach  to  reduce  the  acquisition  time  while  maintaining  the 
4DCBCT  image  quality,  or  from  another  point  of  view,  to 
enhance  the  current  4DCBCT  image  quality  while  using 
short  acquisition  time.  The  approach  developed  in  this  article 
is  essentially  to  reduce  the  view-aliasing  artifacts  resulted 
from  insufficient  sampling.  To  increase  the  sampling  rate  in 
each  individual  phase,  the  projections  of  all  other  phases  can 
be  incorporated  into  the  reconstruction  process.  To  do  so,  a 
motion  model  linking  the  data  of  different  phases  is  required, 
which  can  be  derived  from  deformable  registrations  of  the 
4DCBCT  phases.  Two  other  registration  strategies  of  intro¬ 
ducing  an  additional  diagnostic  CT  image  and  mapping  be¬ 
tween  image  space  and  projection  space  are  also  investi¬ 
gated.  The  proposed  approach  is  evaluated  with  numerical 
simulations.  Since  the  view-aliasing  artifacts  are  prominent 
only  in  the  axial  view  and  for  simplicity,  we  start  with  fan- 
beam  (2D)  geometry  to  study  the  performance  of  the  ap¬ 
proach,  a  clinical  case  of  a  lung  cancer  patient  is  then  pre¬ 
sented. 
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Fig.  1.  Flow  chart  of  the  proposed  reconstruction  process  for  enhanced 
4DCBCT  imaging. 


II.  METHODS  AND  MATERIALS 

II. A.  Modified  Feldkamp  reconstruction  for  reduced 
view-aliasing  artifacts 

As  described  earlier,  to  form  the  4DCBCT  images,  the 
acquired  projection  data  are  retrospectively  sorted  into  sev¬ 
eral  phase  bins  before  reconstructions.  When  the  acquisition 
time  is  limited,  the  reduced  number  of  projections  per  recon¬ 
struction  will  result  in  undersampling  and  degrade  the  image 
with  view-aliasing  artifacts.  In  order  to  improve  the  sampling 
rate  in  each  phase,  we  need  to  borrow  projection  information 
from  other  phases.  Obviously,  direct  inclusion  of  projections 
of  other  phases  will  cause  the  inconsistency  in  the  data  due 
to  patient  motion.  However,  as  we  have  shown  previously,  if 
the  motion  can  be  accurately  modeled,  a  motion-corrected 
CBCT  reconstruction  can  be  obtained.19  This  reconstruction 
method  is  a  modified  Feldkamp  algorithm,  which  performs 
the  backprojections  of  the  filtered  data  along  deformed  paths 
according  to  the  corresponding  motion  model.  Equivalently 
and  computationally  more  efficiently,  one  can  reconstruct 
each  individual  phase  first,  then  superimpose  all  the  phases 
after  deforming  each  of  them  into  the  same  phase  based  on 
the  given  motion  model.  This  method  can  be  described  in  the 
following  flow  chart  (Fig.  1). 

II. B.  Motion  model  estimation 

To  derive  the  motion  model  (i.e.,  the  points  of  correspon¬ 
dence  between  any  two  phases)  to  be  used  in  the  enhanced 
4DCBCT  reconstruction,  we  have  investigated  the  following 
three  registration  strategies  with  and  without  the  aid  of  an 
additional  3D  CT  image. 

II. B.  1.  Register  artifacts-contaminated  4DCBCT 
images 

The  most  straightforward  way  to  find  the  motion  model  is 
to  register  the  regular  4DCBCT  phases  using  a  deformable 


model.  A  unique  problem  here  is  that  the  images  to  be  reg¬ 
istered  contain  serious  steak  artifacts,  which  may  adversely 
affect  the  accuracy  of  the  derived  motion  model,  and  the 
error  can  propagate  into  the  final  4DCBCT  images.  Never¬ 
theless,  the  process  may  prove  to  be  favorable  because  of 
greatly  reduced  view-aliasing  artifacts.  The  performance  of 
this  method  is  shown  later  with  a  simulation  study  where  the 
registration  was  done  based  on  a  free-form  Spline  (BSpline) 
model.  The  simplicity  and  yet  accuracy  of  the  BSpline 
method  make  it  a  preferred  tool  for  many  clinical 
applications.  "  More  details  about  this  registration  algo¬ 
rithm  can  be  found  in  Ref.  22. 


II.B.2.  Register  planning  CT  to  4DCBCT  phased 
image 

Note  that  in  radiation  therapy,  it  is  often  the  case  that  a 
patient  has  already  had  a  3D  or  4D  CT  scan  for  the  treatment 
planning  before  any  radiation  delivery.  In  such  case,  an 
artifact-free  volumetric  image  can  be  assumed  available  prior 
to  the  4DCBCT  scan  on  the  treatment  day.  Considering  res¬ 
piration  motion,  this  CT  image  should  be  one  phase  of  4D 
CT  images  or  a  breath-hold  CT  scan.  The  idea  is  to  register 
each  phase  of  the  4DCBCT  images  to  this  artifact-free  3D 
image  to  derive  corresponding  deformation  with  respect  to 
the  specific  phase  of  the  CT  image.  The  deformation  of  the 
4DCBCT  phases  with  respect  to  each  other  can  be  subse¬ 
quently  obtained  by  calculating  the  relative  movements. 
Since  only  one  of  the  two  images  being  registered  has  arti¬ 
facts,  it  is  expected  that  the  accuracy  will  be  improved  over 
the  first  procedure.  Again,  the  performance  of  this  method  is 
evaluated  with  the  simulation  study  later  in  this  article. 


II.B.3.  Registration  from  image  space  to  projection 
space 

Although  the  reconstructed  4DCBCT  phase  images  con¬ 
tain  view-aliasing  artifacts,  one  should  realize  that  the  arti¬ 
facts  are  not  originally  present  in  the  raw  data  but  are  gen¬ 
erated  later  during  the  reconstruction  process.  Therefore,  we 
may  avoid  the  drawback  of  using  a  low-quality  reconstructed 
image  in  deriving  the  motion  model  by  registering  the  3D 
planning  CT  directly  to  the  CBCT  raw  data  (after  phase  sort¬ 
ing).  Note  that  the  two  images  being  registered  are  in  differ¬ 
ent  spaces;  in  the  following,  we  present  a  simple  algorithm 
for  such  a  registration  task  using  a  BSpline  model. 

Similar  to  the  conventional  image  registration  using  de¬ 
formable  models,  the  motion  of  the  object  is  described  by 
g(x;t)=f(x+n(x',t)),  where  fix)  represents  the  3D  image 
obtained  from  the  planning  CT  scan,  g(x;t)  is  the  same  ob¬ 
ject  at  a  particular  phase  t  during  the  4DCBCT  scan,  and 
r/(jc;f)  defines  the  corresponding  deformation  for  point  x  at 
that  phase.  With  the  BSpline  model,  the  deformation  is  de¬ 
fined  on  a  sparse,  regular  grid  of  control  points  \n  placed 
over  the  CT  image.  When  a  point  is  not  on  the  grid,  the 
associated  deformation  is  calculated  by  the  BSpline  interpo¬ 
lation  as  follows: 
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u(x;t)  =  2  0n  ,f3(i) 

n 


x-K\ 

a,  r 


(i) 


where  /3f3>(x)  =  /3i3\x) /?<3)(y) /3{3\z)  is  a  separable  cubic 
BSpline  convolution  kernel,  0n  t  are  the  values  for  the  control 
points  at  phase  t,  and  At  controls  the  grid  spacing  in  each 
dimension. 

To  register  the  CBCT  projections  with  the  planning  CT 
image,  we  adopt  the  metric  of  sum  of  squared  difference, 

$(0)  =  2  ^  -  2  Aijg^  =  2  (y(  -  2  Aijfpc  +  u)  j", 


(2) 


where  Y ,  is  the  total  attenuation  along  the  z'th  projection  ray, 
j  is  the  index  of  the  image  voxel,  and  A  is  the  cone-beam 
projection  matrix.  For  simplicity,  the  phase  index  t  is  ne¬ 
glected  here.  Given  the  projection  data  Y  and  planning  CT 
image  /,  the  deformation  field  u  at  a  particular  phase  can  be 
found  by  minimizing  the  above  objective  function.  General- 
purpose  algorithms  can  be  applied  to  solve  the  optimization 
problem,  and  in  this  work  the  LBFSG  algorithm  is  used,”6 
which  requires  the  calculation  of  the  first  derivative  of  the 
cost  function  <F.  It  can  be  shown  that  Eq.  (2)  has  a  closed- 
form  derivative  that  can  be  calculated  explicitly.  For  a  simi¬ 
lar  derivation,  readers  are  referred  to  the  original  work  by 
Zeng  et  al.,  ’  of  which  the  method  B3  can  be  understood 
as  a  variation. 


tion  is  possible,  which  may  affect  the  accuracy  of  phase  bin¬ 
ning  in  4DCBCT.  However,  this  was  not  considered  in  our 
simulations  because  the  focus  of  this  paper  is  the  reduction 
of  the  view-aliasing  artifacts. 

II. D.  Patient  study 

A  lung  cancer  patient  had  a  multiple-gantry-rotation 
(MGR)  4DCBCT  scan  in  our  clinic,  and  his  4D  images  were 
used  to  evaluate  our  proposed  approach.  Due  to  the  limit 
access  of  the  available  data,  only  the  method  described  in 
Sec.  2.1.1  was  tested.  The  4DCBCT  was  acquired  with  a 
Trilogy  system  (Varian  Medical  Systems,  Palo  Alto,  CA). 
For  the  MGR  4D  scan,  the  x-ray  tube  current  was  set  to 
32  mA,  and  four  gantry  rotations  were  performed  at  a  speed 
of  6°/s.  Each  rotation  (slightly  over  360°)  consisted  of  over 
680  projections,  and  the  effective  area  of  each  acquired  pro¬ 
jection  image  was  397.312  mm  X  297.984  mm,  containing 
1024  X  768  pixels.  The  average  breathing  cycle  of  the  patient 
was  4.2  s.  The  4DCBCT  images  were  obtained  by  retrospec¬ 
tive  sorting  of  the  MGR  data  into  6  phases  and  subsequently 
reconstructing  the  rebinned  data.  More  details  about  this 
4DCBCT  technique  can  be  found  in  Ref.  16.  To  evaluate  the 
proposed  4D  image  enhancing  method,  the  4D  phases  were 
registered  to  the  0%  phase  with  a  BSpline  deformable  model 
and  then  superimposed.  The  resultant  phase-0%  image  was 
then  compared  with  the  original  4DCBCT  for  several  axial 
slices  to  demonstrate  the  differences. 


II. C.  Simulation  studies 

The  above  model-based  recontruction  method  is  validated 
with  numerical  simulations  under  2D  fan-beam  geometry 
and  3D  cone-beam  geometry.  The  parameters  for  the  simu¬ 
lations  are  similar  to  the  real  setup  of  a  Varian  Trilogy  ma¬ 
chine  (Varian  Medical  Systems,  Palo  Alto,  CA),  where  the 
focal  length  of  the  fan  beam  or  cone  beam  is  150  cm,  and  the 
radius  of  the  detector  rotation  is  50  cm.  The  fan-beam  simu¬ 
lation  has  1024  samples  spaced  by  0.388  mm,  and  2048  pro¬ 
jection  views  were  simulated  with  projection  angles  evenly 
spanned  over  360°,  while  the  cone-beam  simulation  consists 
of  1024X768  pixels  in  each  of  the  total  1200  projections 
with  a  pixel  size  of  0.388  mm.  The  object  used  in  the  fan- 
beam  simulations  came  from  one  slice  of  the  CT  images  of  a 
patient,  which  contained  512  X  512  pixels  with  the  pixel  size 
of  1.0  mm.  For  the  cone-beam  simulation,  the  80  slices  are 
used  with  slicethickness  of  2.5  mm. 

To  include  the  effect  of  the  respiratory-induced  motion, 
we  imposed  an  artificial  perodically  changing  deformation 
field  on  the  image  during  the  simulations.  Specifically,  a 
maximum  deformation  is  determined  by  registering  the  in¬ 
hale  and  exhale  phases  of  a  patient  4DCT  scan.  The  defor¬ 
mation  of  an  intermediate  phase  is  then  obtained  by  linearly 
scaling  down  the  maximum  deformation  field  with  a  scaling 
factor  s = cos2  {irt/T),  where  t  is  the  acquisition  time,  and  T  is 
the  respiration  period.  For  the  fan-beam  case,  the  motion  is 
restricted  in  the  axial  plane,  while  for  the  cone-beam  geom¬ 
etry,  the  motion  is  fully  3D.  The  period  of  the  simulated 
respiratory  motion  is  4  s.  In  reality,  irregular  patient  respira¬ 


III.  RESULTS 

III.A.  Fan-beam  simulation  study 

The  fan-beam  simulated  projection  data  are  shown  in 
Figs.  2(a)  and  2(b),  where  the  left  is  from  the  static  CT  image 
and  the  right  is  from  the  same  object  with  the  artificial  in¬ 
plane  deformation.  It  is  seen  that  the  object  motion  during 
the  scan  generated  a  large  amount  of  inconsistency  in  the 
projection  data  [Fig.  2(b)].  To  show  its  effect  on  the  resultant 
image,  in  Figs.  2(c)  and  2(d),  the  corresponding  recon¬ 
structed  images  are  compared,  where  artifacts  such  as  blur¬ 
ring,  doubling,  distortion,  and  streak  artifacts  are  observed  in 
the  motion  contaminated  image  [Fig.  2(d)]. 

The  simulated  projections  of  the  continuously  deformed 
image  were  retrospectively  sorted  into  eight  phases.  Five  of 
the  reconstructed  phases  from  the  peak  inspiration  to  the 
peak  expiration  are  illustrated  in  Fig.  3(a),  and  a  zoomed-in 
image  of  phase  1  is  shown  in  Fig.  3(b).  Compared  with  the 
reconstructed  images  in  Figs.  2(c)  and  2(d),  it  can  be  found 
that  the  motion  blurring  is  greatly  reduced,  and  the  boundary 
of  the  primary  tumor  became  much  clearer  after  the  phase 
sorting.  However,  some  low  contrast  structures  (for  example 
the  aorta,  muscle,  etc.)  are  lost  due  to  the  view-aliasing 
streak  artifacts. 

Three  enhanced  images  using  the  proposed  approach  are 
shown  in  Fig.  4  corresponding  to  the  three  motion  model 
derivation  techniques.  The  top  row  in  Fig.  4  shows  the  de¬ 
rived  deformation  between  phase  1  and  phase  4,  where  from 
the  left  to  the  right  are  using  “4D  phase-to-phase  registra- 
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Fig.  2.  Simulated  projections  and  their  corresponding  reconstructed  images.  The  left  column  is  from  the  static  object,  and  the  right  column  is  from  the  same 
object  with  artificial  deformation. 


tion”  (method  Bl),  “CT  image  to  4D  phases”  (method  B2), 
and  “CT  image  to  CBCT  projection”  (method  B3),  respec¬ 
tively.  The  corresponding  reconstructed  images  (at  phase  1) 
are  shown  in  the  bottom  row.  We  find  no  significant  artifacts 
in  the  reconstructed  image,  even  with  certain  errors  in  the 
derived  motion  model  [see  Fig.  4(a)  the  circular  pattern  in 
the  deformation  field  due  to  registering  CBCT  images  con¬ 
taminated  with  view-aliasing  artifacts].  The  image  enhancing 


approach  using  any  of  the  three  motion  models  improved  the 
image  quality  over  the  regular  phase-resolved  image  (Fig.  3), 
and  the  method  as  described  in  B3  by  matching  CT  image 
with  CBCT  raw  data  generated  the  best  final  images  (less 
streak  artifacts  near  the  boundary  of  the  field  of  view),  be¬ 
cause  a  relatively  more  accurate  motion  model  is  applied. 

In  Figs.  5(a)  and  5(b),  we  quantitatively  compared  the 
CBCT  images  by  plotting  vertical  profiles  passing  through 


(b) 

Fig.  3.  (a)  Five  examples  of  the  reconstructed  phases  after  the  phase  sorting  of  the  simulated  projection  data  of  the  moving  object;  (b)  a  zoomed-in  image  of 
the  first  reconstructed  phase  where  it  is  found  that  the  motion  blurring  is  reduced  by  phase  binning,  but  some  low-contrast  structures  are  lost  due  to 
view-aliasing  artifacts. 
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Fig.  4.  Top  row,  from  the  left  to  the  right,  illustrates  the  derived  deformation  field  obtained  by  using  methods  Bl,  B2,  and  B3,  respectively,  and  the  bottom 
row  shows  the  corresponding  reconstructed  images  for  phase  1  using  the  proposed  image  enhancement  technique. 


the  tumor.  It  is  seen  from  Fig.  5(a)  that  the  enhanced 
4DCBCT  (ED_4DCBCT)  images  resulted  from  three  mo- 


Relative  Position  (mm) 

<b) 

Fig.  5.  Comparison  of  vertical  profiles  passing  through  the  tumor. 


tion  models  have  similar  accuracy  in  terms  of  the  CT  num¬ 
bers.  The  difference  from  the  phantom  (ground  truth)  is 
mainly  due  to  the  filtered  backprojection  reconstruction  pro¬ 
cess.  In  Fig.  5(b),  the  enhanced  4DCBCT  (using  motion 
model  Bl)  is  compared  with  3DCBCT  and  regular  4DCBCT. 
It  can  be  seen  that  the  ED_4DCBCT  results  in  less  noise 
(compared  with  regular  4DCBCT)  while  maintaining  the 
correct  tumor  profile. 

Furthermore,  the  contrast-to-noise  ratio  (CNR)  as  an  im¬ 
age  metric  is  compared  for  the  CBCT  images,  defined  as 
CNR=|(5>-(50)|  /  cr,  where  (S)  and  (So)  denote  the  average 
signal  within  a  region  of  interest  for  the  tumor  and  soft  tis¬ 
sue,  respectively,  and  cr  denotes  the  standard  deviation 
within  the  tumor.  The  calculated  CNR  for  regular  4DCBCT 
is  1.186;  for  enhanced  4DCBCT,  the  CNRs  are  5.487,  5.584, 
and  5.589  for  motion  models  Bl,  B2,  and  B3,  respectively. 

III.B.  Cone-beam  simulation  study 

Figure  6  shows  the  4D  cone-beam  simulation  results  at 
respiratory  phase  1 .  From  the  top  to  the  bottom  rows  are  the 
phantom  images,  the  motion  blurred  3DCBCT  images,  the 
regular  4DCBCT  images,  and  the  enhanced  4DCBCT  im¬ 
ages,  respectively.  From  left  column  to  the  right  column  are 
the  axial,  coronal,  and  sagittal  views.  The  motion  model  ap¬ 
plied  to  the  enhanced  4DCBCT  image  reconstruction  was 
derived  using  method  B2  by  registering  each  phase  of  the 
regular  4DCBCT  to  the  peak-inhale  phantom  images.  Again, 
dramatic  improvement  [Fig.  6(d)  over  Fig.  6(c)]  is  observed 
by  using  the  proposed  approach  for  the  cone-beam  study. 
The  CNRs  for  the  regular  4DCBCT  and  the  enhanced 
4DCBCT  are  0.92  and  4.43,  respectively. 

III.C.  Patient  study 

The  enhanced  4DCBCT  images  were  obtained  for  the  pa¬ 
tient  using  the  proposed  approach  and  the  motion  model  de- 
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Fig.  6.  Simulation  study  using  cone-beam  geometry.  From  left  to  right  are  the  images  of  axial,  coronal,  and  sagittal  views,  respectively.  Row  (a)  shows  the 
4D  phantom;  row  (b)  is  the  3DCBCT  showing  the  motion  blurring  artifacts;  row  (c)  is  the  regular  4DCBCT  images  showing  the  view-aliasing  artifacts;  and 
row  (d)  is  the  enhanced  4DCBCT  images. 


rived  with  method  Bl.  In  Fig.  7  the  resultant  image  at  phase 
1  is  shown  along  with  the  image  of  the  same  phase  before 
the  enhancement.  Because  MGR  technique  was  used  for  the 
4D  CBCT  scan  for  the  patient,  the  view-aliasing  artifacts  in 
the  original  4DCBCT  images  are  not  prominent.  The  im¬ 
provement  by  using  the  proposed  approach  is  not  dramatic  as 
the  previous  simulation  studies.  The  CNRs  are  3.82  and  3.97 
for  the  regular  4DCBCT  and  enhanced  4DCBCT,  respec¬ 
tively.  It  is  found  that  the  proposed  technique  results  in  a 
smoother  image  in  the  uniform  region  leading  to  a  higher 
CNR.  However,  it  is  also  noticed  that  the  process  may  reduce 
the  spatial  resolution.  The  trade-off  between  the  CNR  and 
the  spatial  resolution  using  the  proposed  model-based  recon¬ 
struction  technique  requires  further  systematic  study,  which 
is  beyond  the  scope  of  this  work. 

IV.  DISCUSSION  AND  CONCLUSION 

One  goal  of  4DCBCT  is  to  remove  the  respiratory  motion 
induced  artifacts  in  the  reconstructed  images  to  increase  the 
accuracy  of  target  localization  and  dose  calculation.  Differ¬ 
ent  from  conventional  4DCT, 29-31  4DCBCT  acquisition  with 
an  onboard  imager  usually  has  a  much  slower  speed  of  about 
1  min  per  rotation,  and  due  to  the  limited  total  scan  time  in 
practice,  undersampling  is  often  an  issue  resulting  in  se¬ 
verely  degraded  4D  images,  which  sometimes  are  even 
worse  than  the  motion  contaminated  3D  images.  The  ap¬ 
proach  proposed  in  this  article  reconstructs  the  image  of  any 
particular  phase  by  introducing  additional  information  from 


other  phases,  effectively  increasing  the  sampling  rate  and 
improving  the  image  quality.  In  general,  the  approach  relies 
on  the  accuracy  of  the  motion  model  that  relates  different 
projections.  We  have  demonstrated  with  2D  simulations  that 
the  best  way  to  derive  the  motion  model  is  to  register  the 
artifact-free  CT  images  with  the  CBCT  projections,  if  a  plan¬ 
ning  CT  scan  has  been  performed  prior  to  the  CBCT  scan. 
Though  the  registration  accuracy  in  3D  geometry  may  be 
different  from  2D  geometry,  these  basic  principles  and  rela¬ 
tive  performances  observed  in  this  article  are  expected  true. 
As  an  example,  the  real  patient  case  shown  in  this  article 
demonstrated  the  feasibility  of  the  proposed  method  under 
the  worst  condition  (low  view-aliasing  artifacts  in  the  4D 
images  leave  little  space  for  further  improvement,  and  low 
accuracy  in  motion  derivation  by  using  method  B 1  limits  the 
performance  of  the  proposed  approach). 

As  it  is  known,  the  radiation  dose  has  been  one  of  the 
major  concerns  in  CBCT  imaging  because  of  the  repeated 
use  for  a  given  patient.  With  the  proposed  approach,  it  is 
possible  to  lower  the  x-ray  tube  current,  hence  the  radiation 
exposure  for  the  4DCBCT  scan  while  maintaining  a  mean¬ 
ingful  4D  image,  because  introducing  additional  information 
of  other  phases  during  reconstruction  effectively  increases 
the  photon  statistics.  This  is  similar  to  the  case  of  4D  fan- 
beam  CT  studied  earlier  by  our  group.  The  efficacy  of  the 
approach  in  improving  the  trade-off  between  the  controlled 
radiation  exposure  and  the  resultant  CBCT  image  quality  is 
being  investigated. 
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Fig.  7.  4DCBCT  images  of  a  patient.  The  left  column  shows  three  arbitrary  slices  of  the  peak  inspiration  phase,  and  the  right  column  is  the  corresponding 
images  after  the  proposed  enhancement  with  motion  model  derived  by  using  method  B 1 . 


It  should  be  noted  that  the  proposed  4DCBCT  enhancing 
technique  needs  a  computation  of  the  deformation  field  be¬ 
fore  the  image  reconstruction,  which  takes  extra  time  com¬ 
pared  with  regular  CBCT  imaging.  Currently,  a  3D  image- 
to-image  registration  with  an  image  size  of  512X512X64 
will  take  approximately  30  min.  The  image-to-projection 
registration  method  generally  takes  longer  because  of  the  cal¬ 
culation  of  the  forward  projections.  The  computational  time 
makes  it  difficult  for  online  applications  in  the  treatment 
room.  However,  the  development  of  techniques  in  speeding 
up  the  algorithms  as  well  as  the  computer  hardware  is  highly 
promising,  which  we  believe  will  make  the  4DCBCT  imag¬ 
ing  fast  and  reliable  in  the  near  future.  In  addition,  the 

4DCBCT  image  quality  may  be  further  improved  by  the  ad- 

33—35 

vanced  reconstruction  algorithm  and  helical  scan  mode." 

In  conclusion,  we  have  developed  a  novel  approach  to 
enhance  the  4DCBCT  image  quality.  The  enhancement  is 
achieved  by  increasing  the  angular  sampling  rate  of  each 
phase  of  the  4DCBCT  data  using  information  from  other 
phases,  thus  the  view-aliasing  artifacts  commonly  seen  in  the 
limited-time  scan  4DCBCT  images  are  eliminated  or  signifi¬ 
cantly  reduced.  The  approach  opens  new  opportunities  for 


4DCBCT  to  be  a  more  useful  tool  that  may  substantially 
improve  the  current  cancer  management  in  radiation  oncol¬ 
ogy- 
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The  purpose  of  this  work  is  to  develop  a  novel  strategy  to  automatically  map  organ  contours  from 
one  phase  of  respiration  to  all  other  phases  on  a  four-dimensional  computed  tomography  (4D  CT). 

A  region  of  interest  (ROI)  was  manually  delineated  by  a  physician  on  one  phase  specific  image  set 
of  a  4D  CT.  A  number  of  cubic  control  volumes  of  the  size  of  ~  1  cm  were  automatically  placed 
along  the  contours.  The  control  volumes  were  then  collectively  mapped  to  the  next  phase  using  a 
rigid  transformation.  To  accommodate  organ  deformation,  a  model-based  adaptation  of  the  control 
volume  positions  was  followed  after  the  rigid  mapping  procedure.  This  further  adjustment  of 
control  volume  positions  was  performed  by  minimizing  an  energy  function  which  balances  the 
tendency  for  the  control  volumes  to  move  to  their  correspondences  with  the  desire  to  maintain 
similar  image  features  and  shape  integrity  of  the  contour.  The  mapped  ROI  surface  was  then 
constructed  based  on  the  central  positions  of  the  control  volumes  using  a  triangulated  surface 
construction  technique.  The  proposed  technique  was  assessed  using  a  digital  phantom  and  4D  CT 
images  of  three  lung  patients.  Our  digital  phantom  study  data  indicated  that  a  spatial  accuracy  better 
than  2.5  mm  is  achievable  using  the  proposed  technique.  The  patient  study  showed  a  similar  level 
of  accuracy.  In  addition,  the  computational  speed  of  our  algorithm  was  significantly  improved  as 
compared  with  a  conventional  deformable  registration-based  contour  mapping  technique.  The  ro¬ 
bustness  and  accuracy  of  this  approach  make  it  a  valuable  tool  for  the  efficient  use  of  the  available 
spatial-tempo  information  for  4D  simulation  and  treatment.  ©  2007  American  Association  of 
Physicists  in  Medicine.  [DOI:  10.1118/1.2780105] 

Key  words:  CT,  image  registration,  contour  mapping,  IGRT 


I.  INTRODUCTION 

A  longstanding  question  in  radiation  therapy  is  how  to  accu¬ 
rately  and  efficiently  segment  a  region  of  interest  (ROI)  such 
as  a  tumor  target  volume  or  a  critical  structure. I_f>  In  spite  of 
intense  research  efforts  in  the  past  few  decades,  ROI  seg¬ 
mentation  remains  a  time  consuming  task  in  treatment  plan¬ 
ning.  In  most  cases,  the  segmentation  is  performed  manually 
in  a  slice-by-slice  fashion,  creating  a  strong  need  for  auto¬ 
mated  segmentation  tools  in  the  clinics.  The  introduction  of 
four-dimensional  computed  tomography  (4D  CT)  in  radia¬ 
tion  oncology  practice  further  amplifies  this  need  as  the  num- 
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ber  of  images  to  be  segmented  is  increased  dramatically. 
Generally,  a  4D  CT  scan  consists  of  5-10  sets  of  three- 
dimensional  (3D)  CT  images,  each  representing  the  patient 
anatomy  at  a  specific  phase  of  respiration.  For  4D  radiation 
therapy  applications,  it  is  labor  intensive  to  follow  the  3D 
approach  of  manual  segmentation  due  to  the  immense  work¬ 
load  associated  with  this  process. 

A  natural  way  to  deal  with  the  4D  segmentation  problem 
is  to  start  with  a  known  set  of  contours  for  a  selected  phase 
and  map  these  contours  onto  all  other  phases.  ROIs  for  the 
selected  phase  are  first  manually  contoured,  similar  to  that 
done  in  treatment  planning  based  on  3D  CT  image  data  sets. 
The  mapping  procedure  can  be  accomplished  with  the  aid  of 
a  computer  algorithm  that  registers  an  arbitrary  point  on  the 
selected  phase  to  the  corresponding  points  on  all  other 


phases.  While  conceptually  simple,  the  implementation  of 
this  idea  is  not  straightforward.  An  intelligent  algorithm  ca¬ 
pable  of  providing  accurate  point-to-point  correspondence 
between  the  phased  images,  or  at  least  between  points  within 
the  ROIs,  is  the  key  to  the  success  of  this  approach.  Various 
studies  have  investigated  algorithms  to  automatically  map 
contours  using  deformable  image  registration  and  surface 
mapping  techniques.  One  method  is  based  on  a  deformable 
registration  model. 16-1 8  This  method  has  limited  accuracy, 
especially  in  the  regions  proximate  to  the  interfaces  of  dif¬ 
ferent  organs,  and  is  brute-force  in  nature,  which  entails  a 
large  amount  of  computations.  In  reality,  contour  mapping  is 
a  regional  problem  and  a  global  association  of  the  phase- 
based  images  is  neither  necessary  nor  efficient.  Surface  map¬ 
ping  achieves  the  stated  goal  of  contour  transformation  by 
iteratively  deforming  the  ROI  contour-extended  surface  until 
the  optimal  match  with  the  reference  is  reached.”’  Nu¬ 
merous  surface  mapping  techniques,  such  as  spatial  parti¬ 
tioning,  principal  component  analysis,  conformal  mapping, 
rigid  affine  transformation,  deformable  contours,  and  warp¬ 
ing  based  on  the  thin-plate  spline  (TPS),  have  been  devel¬ 
oped  over  the  years  and  the  end  point  of  all  of  these  tech¬ 
niques  is  a  mapping  between  topological  components  of  the 
input  surfaces  that  allow  for  transfer  of  annotations.  This 
type  of  computation  is  inherently  more  efficient  in  compari¬ 
son  with  the  deformable  model-based  approaches,  but  it  suf- 
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fers  from  the  fact  that  the  resultant  mapping  heavily  depends 
on  the  model  used  and  the  fact  that  the  model  parameters  in 
the  calculations  are  not  physically  transparent. 

In  this  work,  we  combine  the  useful  features  of  the  two 
different  types  of  techniques  for  contour  mapping.  Our  work¬ 
ing  hypothesis  is  that  information  contained  in  the  boundary 
region  is  often  sufficient  to  guide  the  contour  mapping  pro¬ 
cess  without  relying  on  the  use  of  an  ad  hoc  surface  deform¬ 
ing  model.  In  the  proposed  technique,  the  neighborhood  in¬ 
formation  of  the  contour  surface  is  captured  by  a  series  of 
small  cubic  (or  other  shaped)  control  volumes  placed  around 
the  surface.  The  collection  of  the  centers  of  the  control 
volumes  is  a  representation  of  the  contour  surface.  The  ROI 
contour  mapping  proceeds  iteratively  under  the  guidance  of 
the  information  contained  in  the  control  volumes.  The  pro¬ 
posed  method  is  illustrated  by  a  digital  phantom  experiment 
and  three  clinical  lung  case  studies.  The  results  suggest  that 
the  technique  is  capable  of  automatically  mapping  contours 
among  the  4D  CT  phases  with  clinically  acceptable  accuracy. 


Fig.  1 .  A  schematic  drawing  of  the  placement  of  control  volumes  (shown  as 
squares)  on  the  manually  segmented  contour  (depicted  curve).  Each  control 
volume  is  typically  ~1  cm  in  cubic  shape  as  shown  in  the  zoomed  image. 
No  interpolated  volumes  are  displayed. 


II.  MATERIALS  AND  METHOD 
II. A.  Software  platform 

The  Insight  Toolkit24  (ITK)  and  the  Visualization 
Toolkit25  (VTK),  which  are  open  source  cross-platform 
C++  software  toolkits  and  are  freely  available  for  research 
purposes,  were  used  in  this  study.  A  variety  of  methods  have 
been  programmed  into  the  ITK  platform  for  image  registra¬ 
tion  and  segmentation.  ITK  was  used  for  automatic  mapping 
while  VTK  was  mainly  used  for  3D  visualization  and  con¬ 
tour  construction. 

II. B.  Image  acquisition 

The  4D  CT  image  data  sets  for  three  lung  cancer  patients 
were  acquired  with  a  multislice  helical  CT  scanner  (Discov¬ 
ery  ST,  GE  Medical  System,  Milwaukee,  WI).  The  collected 
data  were  sorted  into  ten  phase  bins.10  The  4D  CT  image  sets 
for  all  patient  studies  were  reconstructed  with  a  2.5  mm  slice 
thickness.  The  size  and  pixel  resolution  of  each  CT  slice  was 
512X512  and  0.98X0.98  mm2,  respectively.  The  4D  CT 
images  were  transferred  through  DICOM  to  a  personal  com¬ 
puter  with  a  Pentium  IV  (2.66  GHz)  processor  for  image 
processing.  One  of  the  phases,  which  is  referred  to  as  the 
template  phase,  was  selected  for  manual  segmentation  of  the 
ROIs.  We  refer  to  the  other  phases  in  the  4D  CT  image  set  as 
the  target  phases  with  corresponding  target  contours.  The 
enrolled  patients  in  this  study  were  under  an  Institutional 
Review  Board  approved  protocol. 

II. C.  Placement  of  control  volumes  along  the  ROI 
contour  surface 

The  task  of  mapping  a  contour  is  to  find  its  corresponding 
location  on  the  target  phase  for  an  arbitrary  point  on  the 
contour  drawn  on  the  template  image.  In  general,  the  image 
feature  surrounding  a  contour  point  can  be  used  as  a  signa¬ 
ture  of  the  point  to  aid  the  search  for  its  corresponding  po¬ 
sition  on  the  target  phase.  In  this  study,  the  image  feature  at 


each  point  was  captured  by  introducing  a  cubic  control  vol¬ 
ume  (~  1  cm  in  size)  centered  at  the  point.  The  manual  ROI 
delineation  was  first  performed  on  the  template  phase  using 
the  Varian  Eclipse  TPS  (Varian  Medical  Systems,  Palo  Alto, 
CA).  Afterwards,  the  contours  were  exported  from  the  TPS 
to  a  local  computer  for  contour  propagation.  The  exported 
contours  are  polygons  on  individual  slices,  and  the  vertices 
of  the  polygons  were  used  as  the  locations  of  control  vol¬ 
umes.  This  is  depicted  in  Fig.  1,  where  the  curve  represents 
the  contour  and  the  squares  represent  the  control  volumes. 
The  center  of  each  control  volume  was  set  at  a  contour  point, 
typically  0.5-1  cm  from  the  next  point.  More  control  vol¬ 
umes  are  placed  through  interpolation  if  the  spacing  between 
two  consecutive  control  volumes  is  greater  than  1  cm.  The 
collection  of  these  points  represents  the  ROI  contour  surface. 

The  contour  mapping  was  carried  out  in  three  steps:  (i) 
mapping  the  introduced  control  volumes  collectively  using  a 
rigid  image  registration  algorithm;  (ii)  iteratively  fine-tuning 
the  3D  positions  of  the  control  volumes  to  determine  the 
deformed  ROI  contour  surface;  and  (iii)  reconstructing  the 
new  contours  by  cutting  the  deformed  surface  slice-by-slice 
along  the  transversal,  sagittal,  or  coronal  direction.  This  pro¬ 
cedure  is  shown  in  Fig.  2  and  described  in  detail  below. 

II. D.  Collective  mapping  of  control  volumes 

After  a  series  of  control  volumes  were  placed  along  the 
segmented  contours  on  the  template  phase,  we  mapped  them 
onto  the  target  phase  collectively  (i.e.,  all  the  control  vol¬ 
umes  were  treated  as  an  entity)  using  a  rigid  image  registra¬ 
tion  algorithm."4  A  feature  of  the  rigid  collective  mapping  is 
that  the  relative  distances  and  orientations  of  the  control  vol¬ 
umes  remain  the  same  during  the  course  of  mapping,  result¬ 
ing  in  an  approximate  ROI  contour  surface  in  the  target 
phase  and  providing  a  good  start  for  further  adjustment  of 
the  contour  shape  to  accommodate  the  organ  deformations. 
We  note  that  using  a  rigid  mapping  is  not  a  necessary  step 
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Fig.  2.  Flow  chart  of  the  automatic  contour  mapping  process  for  4D  CT 
images. 


and  other  methods  capable  of  providing  a  reasonable  initial 
estimate  of  the  ROI  surface  may  also  be  used. 

The  collective  mapping  of  the  control  volumes  was  per¬ 
formed  under  the  guidance  of  a  normal  cross  correlation 
(NCC)  metric24  defined  by 
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where  /a(x,)  is  the  intensity  at  a  point  x,  in  a  control  volume 
indexed  by  a  in  the  target  phase  and  Ip( xj  )  represents  the 
intensity  at  a  point  x'  in  a  control  volume  indexed  by  [i  in 
the  template  phase.  In  Eq.  (1),  x  and  x'  are  related  by  a  rigid 
transformation  T,  where  Tx'=x.  We  calculated  the  transfor¬ 
mation  matrix  T,  which  maps  the  ensemble  of  control  vol¬ 
umes  from  the  template  phase  to  the  target  phase.  The  lim¬ 
ited  memory  Broyden-Fletcher-Goldfarb-Shannon 
algorithm  (L-BFGS)  algorithm26-30  was  used  to  optimize  the 
metric  in  Eq.  (1)  with  respect  to  the  transform  parameters. 
The  details  of  this  algorithm  have  been  described 
elsewhere18  31  and  will  not  be  repeated  here.  During  the 
course  of  the  control  volume  mapping,  an  iterative  calcula¬ 
tion  based  on  the  F-BFGS  algorithm  was  carried  out  until  a 
preset  maximum  number  of  iterations  were  met  or  until  the 
NCC  function  stopped  improving. 


II. E.  Fine  tuning  of  the  mapped  control  volumes 

In  the  absence  of  deformation,  the  ROI  contour  of  the 
target  phase  is  obtained  by  connecting  the  centers  of  the 


rigidly  mapped  control  volumes.  In  a  more  general  case 
where  deformation  does  exist,  the  contour  from  the  rigid 
mapping  serves  as  an  initial  estimate  of  the  ROI  contour  in 
the  target  phase.  Further  positional  adjustment  of  the  control 
volumes  is  needed  to  accommodate  the  deformation  of  the 
ROI.  The  final  positions  of  the  control  volumes,  which  define 
the  ROI  surface,  are  determined  by  balancing  the  self-energy 
and  the  interaction  energy,  which  drive  the  control  volumes 
to  their  corresponding  locations  while  maintaining  similar 
shape  integrity  of  the  contour. 

Mathematically,  the  above  adaptation  process  is  modeled 
by  two  “energy  terms.”  The  first  term  is  referred  to  as  “self¬ 
energy”  and  the  second  term  describes  the  “interaction” 
among  the  control  volumes.  Self-energy  tends  to  drive  a  con¬ 
trol  volume  toward  a  position  where  the  neighborhood  envi¬ 
ronment  resembles  itself  most.  For  each  control  volume  in 
the  template  image,  this  process  was  driven  by  the  NCC 
between  the  volume  and  its  corresponding  locations  in  the 
target  image.  The  interaction  term  among  the  control  vol¬ 
umes  intends  to  maintain  the  integrity  of  the  control  volume 
cluster  as  a  whole  and  prevents  any  unrealistic  control  vol¬ 
ume  configuration  from  happening.  The  interaction  energy 
term  is  described  by 

^interaction  —  ^template  —  ^target’  (^) 

where  Mtemplate  is  defined  as  the  correlation  between  a  con¬ 
trol  volume  in  the  template  phase  and  the  neighboring  con¬ 
trol  volumes  in  the  target  phase,  ,'V/target  represents  the  corre¬ 
lation  between  the  control  volume  in  the  target  phase  and  the 
neighboring  volumes  in  the  template  phase.  These  two  terms 
take  into  account  the  neighborhood  environment  of  the  con¬ 
trol  volume  being  adjusted  and  apply  a  constraint  on  the 
possible  form  of  the  control  volume  configuration.  In  a 
sense,  the  interaction  energy  in  Eq.  (2)  exerts  a  restoring 
force  when  the  position  of  a  control  volume  is  varied  with 
respect  to  its  neighbors.  The  interaction  force  is  important  in 
preventing  the  control  volumes  from  moving  to  unrealistic 
positions  simply  driven  by  the  self-energy  and  in  retaining 
the  shape  integrity  of  the  ROI  surface.  For  simplicity,  only 
the  adjacent  control  volumes  were  considered  when  comput¬ 
ing  the  interaction  energies.  The  final  position  of  each  con¬ 
trol  volume  was  determined  by  minimizing  the  sum  of  the 
self-energy  and  interaction  energy  terms. 

A  simple  searching  algorithm  was  implemented  to  find  the 
minimum  of  Eq.  (2)  and  thus  the  final  configuration  of  the 
control  volumes.  The  algorithm  was  realized  using  an  ex¬ 
haustive  search  in  the  defined  local  region  around  the  control 
volume  and  was  based  on  three  dimensional  basis.  For  each 
rigidly  mapped  control  volume  on  the  target  phase,  we  de¬ 
fined  a  small  region  of  ~3  cm  around  the  central  point  in 
the  volume.  The  search  region  of  3  cm  was  chosen  primarily 
to  accommodate  maximum  deformation  of  tissue.  3  cm  is  a 
very  conservative  value  because  the  tissue  deformation  in 
two  adjacent  phases  barely  exceeds  this  value  (the  maximum 
deformation  occurring  between  the  inhale  and  exhale  phases 
may  reach  ~3  cm).  Equation  (2)  was  then  minimized  to 
determine  the  optimal  control  volume  location  within  this 
region.  This  process  was  repeated  for  each  separate  control 


Medical  Physics,  Vol.  34,  No.  10,  October  2007 


4026 


Chao  et  al.:  Automated  contour  mapping  for  4D  radiation  therapy 


4026 


volume.  Due  to  the  fact  that  the  interaction  energy  term  de¬ 
mands  the  information  from  the  neighboring  volumes,  ad¬ 
justing  the  location  of  each  control  volume  changes  the  po¬ 
sition  of  the  previously  fine-tuned  volume.  Therefore,  several 
cycles  of  this  process  were  needed  to  obtain  the  truly  optimal 
positions  of  all  the  control  volumes.  We  found  (see  digital 
phantom  study  in  results  section)  that  2-3  complete  adjust¬ 
ment  cycles  were  adequate  to  find  the  optimal  control  vol¬ 
ume  configuration. 

II. F.  Reconstruction  of  target  contours 

Upon  completion  of  the  control  volume  mapping  and  ad¬ 
aptation,  the  centers  of  each  control  volume  were  identified. 
A  ROI  surface  was  constructed  with  these  central  points  us¬ 
ing  a  triangulated  surface  construction  technique  which  uses 
marching  cubes  method  and  triangular  surface  decimation 
of  VTK.-3"”  The  intersection  of  the  surface  with  each  CT 
slice  was  superimposed  on  top  of  the  image  in  order  to  vi¬ 
sualize  the  contour  in  a  conventional  fashion.  These  contours 
can  be  exported  as  ASCII  or  DICOM-RT  format  for  treat¬ 
ment  planning. 

II. G.  Evaluation  of  the  algorithm  and  case  study 
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Fig.  3.  Validation  process  of  proposed  algorithm  using  a  digital  phantom, 
(a)  Original  image  and  contour;  (b)  Artificially  deformed  image  and  contour 
(ground  truth  contour)  together  with  contour  after  fine  tuning.  These  two 
contours  are  almost  indistinguishable.  The  contour  after  rigid  mapping  is 
also  shown. 


inspection  is  a  cnvenient  way  for  rapid  assessment  of  a  seg¬ 
mentation  calculation,  especially  in  a  case  where  the  ground 
truth  contours  do  not  exist. 

III.  RESULTS 

III. A.  Digital  phantom  study 


Evaluation  of  a  contour  mapping  algorithm  is  a  difficult 
task  because  of  the  general  lack  of  a  gold  standard  for  com¬ 
parison.  The  proposed  control  volume  based  4D  contouring 
technique  was  first  evaluated  with  a  digital  phantom  experi¬ 
ment.  In  this  study,  a  thoracic  CT  image  was  deformed  in¬ 
tentionally  using  a  known  deformation  matrix  which  was 
introduced  by  drifting  the  positions  of  each  control  volume 
along  the  contour  with  known  sizes.  Specifically,  for  a  con¬ 
trol  volume  i,  we  assign  the  following  displacements: 

D(x,y)  =  i  *  0.01  and  D{£)  =  i  *  0.001  if  i 

=  0,1,...,  N/2  (3) 

or 

D(x,y)  =  (N  -  i)  *  0.01  and  D(z)  =  (N- i)  *  0.001  if  i 

=  N/2  +  1 ,  N/2  +  2,  ...  A,  (4) 

where  N  is  the  total  number  of  control  volumes  on  the  ROI 
contour.  The  artificially  deformed  image  serves  as  a  “new 
breathing  phase”  relative  to  the  original  image.  The  “ground 
truth”  lung  surface  in  the  target  phase  was  attainable  by 
transforming  the  original  contours  with  the  same  deforma¬ 
tion  matrix.  The  manually  delineated  lung  contours  in  the 
original  image  were  also  mapped  using  the  proposed  novel 
control  volume  based  approach.  A  comparison  of  the  mapped 
lung  surface  with  the  ground  truth  allowed  us  to  quantita¬ 
tively  assess  the  success  of  the  proposed  approach. 

The  control  volume  based  contour  mapping  technique 
was  also  applied  to  three  4D  CT  patient  scans.  A  physician 
manually  delineated  the  lungs  and  GTV  for  each  scan  on  a 
selected  phase  (for  example,  the  exhale  phase).  These  con¬ 
tours  were  then  mapped  onto  all  other  phases  using  the  pro¬ 
posed  technique.  Visual  inspection  was  used  to  evaluate  the 
three  patient  studies.  Although  it  is  less  quantitative,  visual 


The  thoracic  CT  images  before  and  after  the  intentionally 
introduced  deformation  are  shown  in  Fig.  3.  The  manually 
delineated  lung  contour  is  shown  in  Fig.  3(a),  while  the 
ground  truth  contour  obtained  by  transforming  the  manual 
contour  with  the  known  transformation  matrix  is  plotted  in 
Fig.  3(b).  The  contour  from  the  proposed  approach  is  shown 
in  the  same  figure.  To  be  comprehensive,  the  rigidly  mapped 
contour,  which  served  as  the  start  of  the  model-based  refine¬ 
ment,  is  also  plotted  in  Fig.  3(b).  Clearly,  it  is  significantly 
deviated  from  the  lung  boundary.  Visual  inspection  of  the 
ground  truth  contour  and  contour  after  refinement  in  Fig. 
3(b)  indicated  that  the  gold  standard  and  the  mapped  con¬ 
tours  are  similar,  despite  the  fact  that  the  intentionally  intro¬ 
duced  deformation  field  is  quite  large  (the  displacement  of 
the  some  of  the  voxels  is  as  large  as  2.0  cm).  The  mean  and 
maximum  separations  between  the  two  sets  of  contours  were 
found  to  be  1.5  and  2.5  mm,  respectively. 

To  better  understand  the  mapping  process,  we  examined 
the  behavior  of  each  energy  term  (metrics)  during  the  course 
of  the  contour  mapping.  In  Fig.  4(a)  we  plotted  the  values  of 
self-energy,  interactive  energy,  and  total  energy  of  one  of  the 
control  volumes  in  the  process  of  positional  adjustment  in 
the  surrounding  area  of  the  control  volume  in  the  deformed 
image.  For  this  particular  control  volume,  the  minimum  of 
Eq.  (2)  was  reached  after  searching  the  first  30  points.  For 
control  volumes  located  in  regions  where  the  deformation  is 
large,  more  searching  points  may  be  required  to  reach  the 
minimum.  It  is  also  interesting  to  show  how  the  average 
metric  value  of  all  the  control  volumes  evolved.  In  Fig.  4(b) 
the  metric  value  is  depicted  as  a  function  of  the  calculation 
process.  The  first  point,  corresponding  to  step  0,  is  the  metric 
before  mapping.  Step  1  refers  to  the  system  after  rigid  map¬ 
ping.  The  metric  at  this  point  is  decreased  but  does  not  reach 
the  minimum.  The  third  point  shows  the  metric  after  each  of 
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Fig.  4.  (a)  Energies  (self,  interaction,  and  total  energy)  as  a  function  of  the  searching  points  in  the  surrounding  area  of  a  randomly  selected  control  volume 
in  the  target  image  (see  text  for  details),  (b)  The  metric  value  versus  step  during  the  course  of  the  contour  mapping.  The  metric  values  are  the  average  NCC 
of  all  the  involved  control  volumes. 


the  control  volumes  is  fine  tuned  sequentially.  To  obtain  the 
optimal  contours  two  more  cycles  of  fine  tuning  were  per¬ 
formed;  however,  the  improvement  with  further  iterations 
was  not  significant. 

III. B.  Patient  studies 

The  proposed  technique  was  assessed  using  4D  CT  im¬ 
ages  of  three  lung  cancer  patients.  Figure  5  shows  the  manu¬ 
ally  delineated  contours  on  the  4D  CT  scans  of  two  lung 
cancer  patients  for  lung  and  GTV,  respectively.  The  CT  im¬ 
ages  and  contours  are  displayed  in  axial,  coronal,  and  sagittal 
views.  Contours  were  drawn  on  the  exhale  phases  (phase  1) 
for  patient  1  (top  row  in  Fig.  5).  The  GTV  contour  for  patient 
3  was  delineated  on  the  inhale  phase  (phase  5)  as  shown  in 
Fig.  5  bottom  row.  For  patients  1  and  2  the  mapping  of  the 
lung  contours  was  studied.  For  the  third  patient,  GTV  con¬ 
tour  mapping  was  investigated. 

The  lung  contours  before  and  after  the  control  volume- 
based  mapping  for  the  first  patient  are  presented  in  Fig.  6.  To 


(a)  (b)  (c) 


Fig.  5.  The  template  CT  images  and  manually  segmented  contours  for  two 
lung  cancer  patients  (top  row:  patient  1;  bottom  row:  patient  3)  on  axial, 
coronal,  and  sagittal  views.  Right  lung  of  patient  1  and  GTV  of  patient  3 
were  manually  delineated,  respectively. 


better  visually  evaluate  the  results,  the  CT  images  and  con¬ 
tours  for  this  patient  are  displayed  in  axial,  coronal,  and  sag¬ 
ittal  views  as  in  left,  middle,  and  right  columns,  respectively. 
Target  phases  4,  7,  and  10  are  selected  to  demonstrate  the 
results  as  in  the  top,  middle,  and  bottom  rows  in  Fig.  6.  The 
template  contours  and  rigidly  mapped  contours  were  shown 
and  significantly  deviated  from  the  ROI  boundary.  The  final 
contours  after  the  model-based  adaption  are  presented  as 
well. 

When  the  lung  deformation  is  small,  for  example,  in 
phase  10  of  patient  1  [see  Figs.  6(g)-6(i)],  the  rigidly 
mapped  contours  closely  resemble  the  target  contours.  Fine 
tuning  merely  provides  limited  adjustment  of  the  contours. 
For  phases  with  large  deformations  (e.g.,  phases  of  4  and  7 


Fig.  6.  Axial,  coronal,  and  sagittal  views  of  CT  images  along  with  lung 
contours  for  the  first  patient.  Top  row:  phase  4;  middle  row:  phase  7;  bottom 
row:  phase  10.  Rigidly  mapped  and  target  contours  are  displayed.  Template 
contour  is  also  overlaid  on  displayed  phases. 
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Fig.  7.  Axial,  coronal,  and  sagittal  views  of  CT  images  along  with  GTV 
contours  for  the  third  patient.  Top  row:  phase  2;  middle  row:  phase  4; 
bottom  row:  phase  7.  Rigidly  mapped  and  target  contours  are  displayed. 
Template  contour  is  also  overlaid  on  displayed  phases. 


of  this  patient),  the  model-based  adjustment  after  the  rigid 
contour  mapping  plays  an  important  role  in  obtaining  opti¬ 
mal  contour  arrangement.  In  case  of  the  lung  contours,  our 
analyses  of  the  patient  data  indicate  that  the  accuracy  of  our 
contour  mapping  technique  is  within  one  pixel  in  the  supe¬ 
rior  part  of  the  lung,  where  respiration  induced  deformation 
is  small.  In  other  parts,  the  average  error  on  the  target  con¬ 
tours  is  estimated  to  be  less  than  2.5  mm. 

The  mapped  GTV  contours  for  phases  2,  4,  and  7  in  the 
study  of  patient  3  are  shown  in  Fig.  7.  The  representation  of 
the  contours  is  the  same  as  in  Fig.  6.  For  phases  4  and  7  the 
deformation  was  relatively  small  and  our  contour  mapping 
algorithm  performed  reliably.  For  phase  2,  which  corre¬ 
sponded  to  the  exhale  phases,  significant  deformation  in  the 
ROIs  was  observed.  Despite  this  less  ideal  case,  our  algo¬ 
rithm  still  worked  well  in  these  phases. 

IV.  DISCUSSION 

As  in  3D  radiation  therapy,  delineation  of  ROIs  in  4D  CT 
images  is  a  necessary  step  for  treatment  planning.34  4D  seg¬ 
mentation  is  required  for  constructing  a  4D  patient  model 
and  for  computing  the  accumulated  dose  of  moving  or  de¬ 
formed  organs.  One  way  to  proceed  is  to  draw  all  contours 
on  one  of  the  phases  of  respiration  and  map  or  propagate 
these  contours  onto  the  remaining  phases.  Various  deform¬ 
able  models  have  been  used  for  ROI  contour  mapping  as 
discussed  in  the  introduction.  In  this  work,  we  take  a  re¬ 
gional  approach  based  on  the  mapping  and  adaptation  of  the 
sparsely  sampled  control  volumes.  Our  approach  takes  ad¬ 
vantage  of  the  imaging  features  surrounding  the  ROI  and 
uses  them  as  guidance  in  searching  for  the  optimal  mapped 
contours  while  considering  the  shape  integrity  of  the  ROI 
surface. 


The  mapping  of  a  point  in  one  image  to  another  is  easily 
achievable  if  a  unique  identifier  or  signature  can  be  tagged  to 
the  point.  Here  the  image  feature  contained  in  a  control  vol¬ 
ume  is  employed  as  a  signature  of  the  point  to  facilitate  the 
process  of  finding  its  corresponding  location  in  the  target 
phases.  The  concept  of  control  volume  was  first  introduced 
by  Schreibmann  and  Xing"  and  its  advantages  for  both 
intra-  and  intermodality  image  registration  have  been  dem¬ 
onstrated.  This  study  represents  a  novel  application  of  the 
concept  to  4D  image  segmentation.  After  a  rigid  mapping  of 
the  contours  from  the  template  phase  to  the  target  phase,  a 
model-based  adaptation  is  performed  to  establish  a  reliable 
association  between  the  ROIs  in  two  phase  specific  image 
sets.  This  adaptation  takes  into  account  the  deformation  of  an 
object  and  ensures  the  overall  integrity  of  the  resultant  ROI 
surface.  The  calculation  is  local  in  nature,  which  improves 
both  computational  efficiency  and  convergence  behavior. 

A  quantitative  comparison  of  different  deformable  regis¬ 
tration  algorithms  is  a  difficult  task  because  of  the  multifac¬ 
eted  and  even  subjective  nature  of  the  problem.  A  unbiased 
and  meaningful  comparison  may  entail  the  efforts  from  mul- 

7C 

tiple  institutions.  Thus  we  defer  the  detailed  comparative 
study  to  the  future.  However,  we  wish  to  emphasize  that  the 
chief  advantage  of  the  proposed  technique  is  its  computa¬ 
tional  efficiency.  Because  of  the  regional  nature  of  the  calcu¬ 
lation,  our  general  finding  is  that  the  new  algorithm  is  at  least 
an  order  of  magnitude  faster  than  the  whole  image  based 
approach.  ’  Enormous  saving  in  the  memory  usage  in  the 
proposed  approach  is  also  self-explanatory  and  highly  desir¬ 
able  feature  in  practice. 

A  common  problem  in  image  segmentation  and  contour 
mapping  studies  is  the  lack  of  quantitative  validation.  In  the 
study  of  Lu  et  al.  ,36  for  example,  the  accuracy  of  a  deform¬ 
able  model-based  contour  mapping  technique  was  evaluated 
purely  based  on  visual  inspection.  The  same  approach  was 
employed  in  many  other  previous  investigations.4’6’16,18  In 
our  study,  in  addition  to  the  visual  evaluation,  a  set  of  digital 
phantom  experiments  was  introduced  to  evaluate  the  success 
of  the  proposed  technique.  By  applying  a  prespecified  defor¬ 
mation  matrix  to  the  original  image,  the  ground  truth  of  the 
contour  propagation  is  readily  known.  Therefore,  the  experi¬ 
ments  provide  a  quantitative  test  of  the  proposed  algorithm. 
In  general,  we  found  that  a  spatial  accuracy  better  than  2.5 
mm  is  achievable  using  our  technique. 

V.  CONCLUSION 

The  development  of  4D  radiation  therapy  involves  the  use 
of  a  large  number  of  images  acquired  at  different  times 
and/or  with  different  modalities.  Clinical  implementation  of 
the  new  IGRT  paradigm  is,  to  a  large  extent,  bottlenecked  by 
the  inability  to  accurately  and  efficiently  register  images  and 
segment  the  ROIs.  In  this  work,  a  mathematical  framework 
for  control  volume-based  contour  mapping  has  been  pro¬ 
posed  for  4D  radiation  therapy.  We  demonstrated  that  the 
information  contained  in  the  boundary  region  is  sufficient  to 
guide  the  contour  mapping  process  without  registering  the 
whole  image  or  relying  on  the  use  of  an  ad  hoc  surface 
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deforming  model.  The  results  showed  that  the  control 
volume-based  contour  mapping  algorithm  is  capable  of  ro¬ 
bustly  and  accurately  mapping  contours  from  one  phase  of  a 
4D  CT  to  the  remaining  phases.  Our  technique  decreases  the 
workload  involved  in  4D  CT  ROI  segmentation  and  provides 
a  valuable  tool  for  the  efficient  use  of  available  spatial-tempo 
information  for  4D  simulation  and  treatment. 
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AUTOMATED  CONTOUR  MAPPING  WITH  A  REGIONAL  DEFORMABLE  MODEL 


Ming  Chao,  Ph.D.,  Tianfang  Li,  Ph.D.,  Eduard  Schreibmann,  Ph.D., 

Albert  Koong,  M.D.,  and  Lei  Xing,  Ph.D. 

Department  of  Radiation  Oncology,  Stanford  University  School  of  Medicine,  Stanford,  CA 

Purpose:  To  develop  a  regional  narrow-band  algorithm  to  auto-propagate  the  contour  surface  of  a  region  of  inter¬ 
est  (ROI)  from  one  phase  to  other  phases  of  four-dimensional  computed  tomography  (4D-CT). 

Methods  and  Materials:  The  ROI  contours  were  manually  delineated  on  a  selected  phase  of  4D-CT.  A  narrow  band 
encompassing  the  ROI  boundary  was  created  on  the  image  and  used  as  a  compact  representation  of  the  ROI  sur¬ 
face.  A  BSpline  deformable  registration  was  performed  to  map  the  band  to  other  phases.  A  Mattes  mutual  infor¬ 
mation  was  used  as  the  metric  function,  and  the  limited  memory  Broyden-Fletcher-Goldfarb-Shanno  algorithm 
was  used  to  optimize  the  function.  After  registration  the  deformation  field  was  extracted  and  used  to  transform 
the  manual  contours  to  other  phases.  Bidirectional  contour  mapping  was  introduced  to  evaluate  the  proposed  tech¬ 
nique.  The  new  algorithm  was  tested  on  synthetic  images  and  applied  to  4D-CT  images  of  4  thoracic  patients  and 
a  head-and-neck  Cone-beam  CT  case. 

Results:  Application  of  the  algorithm  to  synthetic  images  and  Cone-beam  CT  images  indicates  that  an  accuracy  of 
1.0  mm  is  achievable  and  that  4D-CT  images  show  a  spatial  accuracy  better  than  1.5  mm  for  ROI  mappings  be¬ 
tween  adjacent  phases,  and  3  mm  in  opposite-phase  mapping.  Compared  with  whole  image-based  calculations, 
the  computation  was  an  order  of  magnitude  more  efficient,  in  addition  to  the  much-reduced  computer  memory 
consumption. 

Conclusions:  A  narrow-band  model  is  an  efficient  way  for  contour  mapping  and  should  find  widespread  applica¬ 
tion  in  future  4D  treatment  planning.  ©  2008  Elsevier  Inc. 

Deformable  model,  Image  registration,  Contour  mapping,  IGRT. 


INTRODUCTION 

Segmentation  of  a  region  of  interest  (ROI),  such  as  a  tumor 
target  volume  or  a  sensitive  structure,  is  an  important  but 
time-consuming  task  in  radiotherapy  (1-6).  With  the  emer¬ 
gence  of  four-dimensional  (4D)  imaging  and  adaptive  radio¬ 
therapy,  the  need  for  efficient  and  robust  segmentation  tools 
is  even  increasing  (7-1 1).  Because  of  dramatically  increased 
numbers  of  images,  it  becomes  impractical  to  manually  seg¬ 
ment  the  ROIs  slice  by  slice  as  in  current  three-dimensional 
radiotherapy  practice.  A  natural  solution  to  the  4D  computed 
tomography  (4D-CT)  segmentation  problem  is  to  delineate 
the  ROIs  on  a  selected  phase  and  then  propagate  the  contours 
onto  other  phases  using  a  mathematical  model.  Along  this 
line,  deformable  model-based  contour  mapping  has  been 
implemented  by  a  few  groups  (12-14).  Although  feasible, 
the  calculation  is  global  in  nature  and  thus  computationally 
intensive.  In  addition,  the  accuracy  of  the  mapped  contours 
may  be  compromised  because  the  registration  may  be 


influenced  unnecessarily  by  the  image  content  distant  from 
the  ROIs,  which  would  otherwise  be  irrelevant  to  the  contour 
mapping  process.  This  is  especially  problematic  when  non¬ 
local  deformable  models,  such  as  thin  plate  spline  and  elastic 
model,  are  used.  In  general,  contour  mapping  is  a  regional 
problem,  and  a  global  association  of  the  phase-based  images 
is  neither  necessary  nor  efficient. 

Surface  mapping  techniques  (15-17)  represent  a  competi¬ 
tive  alternative  to  the  defonnable  model-based  approach. 
The  idea  of  surface  mapping  is  to  obtain  contour  transforma¬ 
tion  by  iteratively  defonning  the  ROI  contour-extended  sur¬ 
face  until  the  optimal  match  with  the  reference  is  found.  The 
calculation  involves  only  the  surface  region  and  is  thus  com¬ 
putationally  efficient.  Numerous  surface  mapping  techniques 
have  been  developed  in  the  past,  which  include,  to  name 
a  few,  spatial  partitioning,  principal  component  analysis, 
conformal  mapping,  rigid  affine  transformation,  deformable 
contours,  and  warping  based  on  the  thin-plate  spline.  All  of 
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these  techniques  are  a  mapping  between  topologic  compo¬ 
nents  of  the  input  surfaces  that  allow  for  transfer  of  annota¬ 
tions.  Although  the  calculations  are  inherently  efficient,  the 
results  depend  heavily  on  the  model  used,  which  may  not 
be  generally  applicable  for  all  clinical  situations  because 
the  ROI  surface  is  multidimensional  and  hardly  modeled 
by  only  a  few  parameters. 

In  this  work,  we  present  a  novel  regional  algorithm  for  ROI 
propagation  among  different  4D-CT  phases.  The  deforma¬ 
tion  of  an  ROI  contour-extended  surface  in  our  algorithm  is 
not  driven  by  an  ad  hoc  surface-based  model  but  instead  by 
the  image  features  in  the  neighborhood  of  the  surface.  The 
underlying  hypothesis  here  is  that  information  contained  in 
the  ROI  boundary  region  is  sufficient  to  guide  the  contour 
mapping  process.  In  the  proposed  algorithm  the  neighbor¬ 
hood  image  features  of  an  ROI  are  captured  by  a  narrow 
band,  which  is  composed  of  all  points  within  two  surfaces 
with  the  signed  distances  of  ±d  from  the  ROI  boundary. 
The  algorithm  is  a  hybrid  of  the  regional  surface-based 
model  and  the  global  deformable  registration-based  ap¬ 
proach.  The  combination  takes  advantage  of  the  desirable 
features  of  each  of  these  two  techniques  and  provides  a  robust 
and  computationally  efficient  contour  propagation  tool  for 
4D  radiotherapy. 

METHODS  AND  MATERIALS 

Software  platform 

The  proposed  contour  mapping  algorithm  was  implemented  using 
the  Insight  Toolkit  (18)  and  the  Visualization  Toolkit  (19),  which 
are  open  source  cross-platform  C++  software  toolkits  sponsored 
by  the  National  Library  of  Medicine. 

Overview  of  the  mapping  process 

Figure  1  depicts  the  overall  contour  mapping  process.  For  a  given 
4D-CT  image  set,  a  selected  phase,  named  the  template  phase,  was 
selected,  and  the  ROIs  were  manually  delineated  by  a  physician.  The 
manually  outlined  contour  was  referred  to  as  the  template  contour.  A 
narrow  band  encompassing  the  template  contour  was  created  (see 
next  section  for  details).  A  deformable  mapping  was  then  carried 
out  to  propagate  the  band  from  the  template  phase  to  other  phases, 
referred  to  as  target  phases.  Upon  successful  mapping  of  the  band, 
the  deformation  field  was  used  to  transform  the  template  contour 
to  the  target  images. 

Narrow-hand  representation  of  ROI  contour 

The  contour  manually  segmented  on  an  axial  slice  of  the  template 
image  has  a  polygon  shape,  and  the  vertices  of  the  polygon  form  the 
basis  for  constructing  the  narrow  band.  As  schematically  shown 
in  Fig.  2,  a  band  with  signed  distances  ±d  was  placed  along  the  tem¬ 
plate  contour.  The  regional  image  features  contained  in  the  band 
function  serve  as  a  “signature”  of  the  contour  and  drive  the  contour 
mapping  process.  The  distance  between  the  neighboring  vertices  on 
the  contour  is  typically  2-10  mm,  depending  on  the  shape  of  the 
contour.  In  generating  the  narrow  band,  we  first  created  cubes 
with  a  side  length  of  2d  around  all  the  vertices,  as  depicted  by  points 
A  and  B  in  Fig.  2.  To  obtain  a  smooth  band,  between  A  and  B  three 
more  cubes,  centered  at  points  C,  D,  and  E,  were  inserted.  Point  C 
was  chosen  to  be  the  middle  point  between  A  and  B,  point  D  the 
middle  between  A  and  C,  and  point  E  the  middle  between  C  and 
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(a) 
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Fig.  1.  Flow  chart  of  narrow  band-based  contour  mapping  proce¬ 
dure.  (a)  Overall  calculation  process,  (b)  Deformable  mapping  pro¬ 
cess  of  the  narrow  band. 

B.  More  interpolated  vertex  points  can  be  introduced  similarly 
when  needed.  Figure  3  illustrates  a  narrow  band  surrounding  the 
lung  boundary  on  the  template  phase  CT  image.  The  light  green 
area  stands  for  the  narrow  band,  and  the  green  curve  is  the  manual 
contour.  The  width  of  the  narrow  band  was  set  to  be  2d  =  15  mm 
in  our  calculations.  To  examine  the  robustness  of  the  proposed  map¬ 
ping  algorithm,  a  variety  of  other  bandwidths,  ranging  from  4  mm 
through  30  mm,  were  also  tested  for  one  of  the  clinical  cases. 

Contour  propagation 

As  illustrated  in  Fig.  1,  the  process  of  contour  mapping  is  essen¬ 
tially  to  waip  the  narrow  band  constructed  above  in  such  a  way 
that  its  best  match  in  the  target  image  is  found.  Mathematically, 
the  mapping  process  of  the  narrow  band  constitutes  an  optimization 
problem,  in  which  a  group  of  transformation  parameters  that  trans¬ 
form  the  points  within  the  band  in  the  template  phase  to  their  homol¬ 
ogous  points  in  the  target  image.  The  warping  of  the  narrow  band  is 
quantified  by  a  metric  function,  which  ranks  a  trial  matching  based  on 
the  “accordance”  level  of  the  image  content  of  the  band  and  its  cor¬ 
respondence  in  the  target  image.  The  calculation  process  is  detailed 
below. 
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Narrow  Band 


Fig.  2.  A  schematic  drawing  of  narrow-band  construction. 


The  input  to  the  contour  mapping  software  includes  the  narrow 
band  and  the  whole  target  image,  which  are  described  by  the  image 
intensity  distributions  Ia(x)  and  Ib(x),  respectively.  It  is  worth  em¬ 
phasizing  that,  even  though  the  whole  target  image  was  used,  only 
fractional  voxels  in  the  target  image  (the  voxels  encompassed  by 
the  band)  are  involved  in  each  iteration  (a  subregion  surrounding 
the  ROI  on  the  target  image  could  be  created  and  used  in  the  calcu¬ 
lation,  but  the  algorithm  converged  so  fast  that  after  two  to  three  it¬ 
erations  the  searching  was  quickly  confined  in  the  neighborhood  of 
the  optimal  solution).  The  narrow  band  acts  as  a  representation  of 
the  ROI  contour.  The  task  is  to  find  the  transformation  matrix, 
T(x),  that  maps  an  arbitrary  point  in  the  band  to  the  corresponding 
point  on  the  target  image  (or  vice  versa)  so  that  the  best  possible  cor¬ 
respondence,  as  measured  by  the  metric  function,  is  achieved.  The 
calculation  proceeds  iteratively.  A  BSpline  deformable  model  is 
used  to  model  the  deformation  of  the  band,  but  other  models  should 
also  be  applicable.  The  spacing  between  the  BSpline  nodes  was  cho¬ 
sen  to  be  approximately  0.5  cm  (smaller  spacing  was  tested,  but  no 
significant  difference  was  found  in  the  final  registration  results).  The 
displacement  of  a  node  ;  is  specified  by  a  vector  X;,  and  the  displace¬ 
ment  vectors  (20)  of  a  collection  of  nodes  characterize  the  tissue 
deformation.  The  displacement  at  a  location  x  on  the  image  is  de¬ 
duced  by  a  BSpline  polynomial  fitting. 

The  Mattes  Mutual  Information  (MMI)  (21)  was  used  as  the  met¬ 
ric  function  for  narrow-band  mapping  (22-25).  The  central  concept 


Fig.  3.  Computed  tomographic  images  with  manual  contours  and  the  narrow  bands  for  patient  1.  The  narrow  bands  are 
shown  in  light  green  and  the  contours  are  green  curves,  (a)  Transverse  view;  (b)  coronal  view;  (c)  sagittal  view. 
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of  mutual  information  (MI)  is  the  calculation  of  entropy.  For  an  im¬ 
age  A,  the  entropy  is  defined  as 


H 


=  -  j  Pa  (a) log  pA(a)da1 


where  pA(a)  (also  called  the  marginal  probability  density  function 
[PDF])  is  the  probability  distribution  of  grey  values  (image  intensi¬ 
ties),  which  is  estimated  by  counting  the  number  of  times  each  grey 
value  occurs  in  the  image  and  dividing  those  numbers  by  the  total 
number  of  occurrences.  Given  two  images,  A  and  B,  their  joint  en¬ 
tropy  is 


H 


J  J  Pab(ci,  b)  log  pAR{a,  b)dadb , 


where  pAB(a,b)  is  the  joint  PDF  defined  by  a  ratio  between  the 
number  of  grey  values  in  the  joint  histogram  (feature  space)  of 
two  images  and  the  total  entries  (26).  The  mutual  information  is  gen¬ 
erally  expressed  as 


MI(A,B)  =  H{A)  +  H(B)  -  H(A,B). 

Mutual  information  measures  the  level  of  information  that  a  ran¬ 
dom  variable  (e.g.,  Ia(x))  can  predict  about  another  random  variable 
(e.g.,  Ib(x)).  Different  from  the  conventional  MI,  whereby  two  sep¬ 
arate  intensity  samples  are  drawn  from  the  image,  the  Mattes  imple¬ 
mentation,  MMI,  uses  only  one  set  of  intensity  to  evaluate  both  the 
marginal  and  joint  PDFs  at  discrete  positions  or  bins  that  uniformly 
spread  within  the  dynamic  range  of  the  images.  Entropy  values  were 
computed  by  summing  over  all  the  bins.  The  number  of  bins  used  to 
compute  the  entropy  in  MMI  metric  evaluation  was  chosen  to  be  30, 
and  the  number  of  spatial  samples  used  was  20,000.  Details  of  MMI 
implementation  can  be  found  in  Mattes  et  al.  (21). 

The  limited  memory  Broyden-Fletcher-Goldfarb-Shanno  algo¬ 
rithm  (L-BFGS)  (27-29)  was  used  to  optimize  the  MMI  metric  func¬ 
tion  with  respect  to  the  displacement  parameters  of  the  nodes,  { x;  j , 
to  find  the  transformation  matrix  T(x)  that  relates  the  points  on 
image  A  and  image  B.  Here  we  just  briefly  show  the  algorithm. 
Starting  from  a  positive  definitive  approximation  of  the  inverse 
Hessian  H0  at  x0,  L-BFGS  derives  the  optimization  variables  by 
iteratively  searching  through  the  solution  space.  At  an  iteration  k, 
the  calculation  proceeds  as  follows:  [I]  determine  the  descent 
direction  pt.  =  — H iiVf(xic)\  [2]  line  search  with  a  step  size 
oik  =  ar*  >'()  (x*  +  apk),  where  a  is  the  step  size  defined  in  the  L- 
BFGS  software  package;  [i]  update  xk+1  =  xk  +  ak  pk  ;  and  [4]  com¬ 
pute  Hk+i  with  the  updated  Hk  . 

At  each  iteration  a  backtracking  line  search  is  used  in  L-BFGS  to 
determine  the  step  size  of  movement  to  reach  the  minimum  of/ 
along  the  ray  xk  +  apk.  For  convergence  a  has  to  be  chosen  such 
that  a  sufficient  decrease  criterion  is  satisfied,  which  depends  on 
the  local  gradient  and  function  value  and  is  specified  in  L-BFGS 
by  the  Wolfe  conditions  (27).  During  the  course  of  optimization, 
the  above  iterative  calculation  based  on  L-BFGS  algorithm  con¬ 
tinues  until  the  following  stopping  criterion  is  fulfilled: 


max(l,  INU 

or  a  pre-set  maximum  number  of  iterations  is  reached.  In  this  study 
we  set  £  =  10®  and  the  iteration  number  to  200,  but  no  more  than  100 
iterations  were  exceeded  in  all  our  calculations  for  the  algorithm  to 
converge. 


Evaluation  of  algorithm  performance 

Evaluation  of  a  contour  mapping  algorithm  is  a  difficult  task  be¬ 
cause  of  the  lack  of  the  ground  truth  for  comparison.  A  straightfor¬ 
ward  means  of  evaluation  is  the  visual  inspection  of  the  mapped 
contours.  In  addition  to  this,  evaluation  based  on  synthetic  images 
(digital  phantoms)  is  also  commonly  used.  The  images  and  existing 
contours  are  distorted  with  preset  deformation  fields.  Because  the 
gold  standard  is  known,  a  direct  comparison  with  the  mapped  con¬ 
tour  is  made  so  as  to  assess  the  propagation  algorithm  quantitatively. 
Beside  these  two  methods,  we  further  performed  a  bidirectional 
mapping  to  evaluate  the  proposed  algorithm.  In  this  test,  the  reverse 
of  the  original  contour  mapping  was  performed:  the  mapped  con¬ 
tours  on  the  target  phase  were  treated  as  the  template  contours  and 
mapped  back  to  the  original  template  phase.  The  contours  so  ob¬ 
tained  were  then  compared  with  the  original  manual  contours,  and 
the  difference  between  the  two  sets  of  contours  was  quantified. 
The  difference  between  the  resultant  and  template  contours  was 
measured  in  terms  of  the  displacements  of  the  vertex  points  on  the 
two  contours.  The  last  yet  pragmatic  evaluation  of  the  algorithm  per¬ 
formance  on  patient’s  study  was  based  on  the  physician’s  manual 
contours. 

Case  study 

Four  thoracic  cancer  patients,  named  as  patient  1, 2,  3,  and  4,  were 
first  used  to  test  the  proposed  algorithm.  These  patients  underwent 
4D-CT  scans.  The  4D-CT  images  were  acquired  with  a  GE  Discov- 
ery-ST  CT  scanner  (GE  Medical  System,  Milwaukee,  WI).  The  col¬ 
lected  data  were  sorted  into  10  phase  bins.  The  ROIs  on  the  template 
phase  were  manually  segmented  by  a  physician.  Specifically,  for  pa¬ 
tients  1  and  2,  the  inhale  phase  was  chosen  for  manual  segmentation, 
and  for  patients  3  and  4,  the  exhale  phase.  Different  ROIs  were  used 
to  better  evaluate  the  algorithm.  Lungs  were  selected  from  patients 
1,  2,  and  3  and  gross  tumor  volume  (GTV)  from  patient  4.  Figure  3 
illustrates  the  manual  contour  and  narrow  band  representation  for 
the  lung  from  patient  1 .  Contour  is  shown  in  the  green  curve  together 
with  the  regional  narrow  bands  (light  green  area)  on  the  transverse, 
coronal,  and  sagittal  views  (Figs.  3a,  3b,  and  3c,  respectively). 

To  further  assess  the  robustness  of  the  proposed  algorithm,  we 
also  carried  out  the  contour  propagation  calculation  from  planning 
CT  to  Cone-beam  CT  (CBCT)  for  a  head-and-neck  case.  The 
CBCT  images  were  acquired  using  the  Varian  Trilogy  system  (Var- 
ian  Medical  Systems,  Palo  Alto,  CA). 

RESULTS 


Convergence  analysis 

To  better  illustrate  the  iterative  process  of  the  contour 
propagation,  in  Fig.  4  the  MMI  metric  as  a  function  of  itera¬ 
tion  step  is  plotted  for  the  narrow  band  mapping  from  the  first 
phase  (inhale  phase)  to  the  other  nine  phases  for  the  first  tho¬ 
racic  patient.  In  all  nine  calculations  it  is  seen  that  the  metric 
value  decreases  monotonically  as  the  iteration  proceeds. 
However,  the  number  of  iterations  needed  for  the  algorithm 
to  find  the  optimal  solution  varies.  It  is  interesting  to  observe 
that,  for  an  “easier”  mapping  whereby  the  deformation  be¬ 
tween  the  two  phases  is  small,  the  number  of  iterations 
required  is  less,  whereas  for  “tougher”  ones  with  larger  dif¬ 
ferences  in  ROI  shapes,  the  required  number  of  iterations  in¬ 
creases  drastically.  Indeed,  from  Fig.  4  it  is  seen  that  the 
minimum  number  of  iterations  required  for  the  metric  to  sat¬ 
urate  occurs  when  mapping  the  phase  1  to  the  adjacent 
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Iteration 

Fig.  4.  Narrow-band  metric  values  as  a  function  of  iteration  step 
when  mapping  the  narrow  band  from  phase  1  to  the  other  nine 
phases  of  the  four-dimensional  computed  tomography. 


phases,  2  and  10.  For  other  mapping,  the  required  iteration  in¬ 
creases  and  reaches  its  largest  value  for  the  “toughest”  map¬ 
ping  between  inhale  and  exhale  (phase  5)  phases. 

In  the  above  analysis,  the  bandwidth  was  set  to  be  15  mm. 
The  performance  of  the  proposed  algorithm  was  also  evalu¬ 
ated  by  varying  the  width  in  the  range  of  4  mm  and  30 
mm.  Specifically,  we  tried  the  widths  of  4  mm,  8  mm,  10 
mm,  15  mm,  20  mm,  and  30  mm.  Our  results  revealed  that, 
when  the  band  was  too  narrow  (e.g.,  4  mm),  the  mapping 
may  fail  locally  at  a  place  not  containing  sufficient  neighbor¬ 
hood  image  features.  The  situation  is  improved  dramatically 
as  the  bandwidth  increases.  For  all  the  clinical  cases  studied 
here,  no  single  failure  was  observed  for  a  width  of  15  mm. 
When  the  width  is  too  large,  the  whole  ROI  will  be  included 
in  the  band.  In  this  situation,  the  mapping  becomes  equiva¬ 
lent  to  registering  the  whole  image  and  the  advantage  of 
the  narrow  band  will  be  overshadowed  by  the  dramatically 
increased  memory  and  computing  costs.  Our  experience  indi¬ 
cates  that  a  width  of  10-15  mm  provides  a  fine  balance 
between  the  computational  accuracy  and  the  associated  cost. 

We  found  that  the  overall  computing  time  was  increased 
by  roughly  an  order  of  magnitude  when  going  from  the  nar¬ 
row  band  approach  to  the  conventional  deformable  model- 
based  contour  mapping,  say,  approximately  3  min  for  narrow 
band-based  mapping  vs.  approximately  25  min  for  whole 
image-based  mapping.  The  dramatically  increased  computer 
memory  requirement  in  the  latter  case  also  posts  a  serious 
problem  when  developing  a  clinically  practical  contour  prop¬ 
agation  method  for  4D  radiotherapy. 


Algorithm  performance  evaluation 

In  addition  to  visual  inspect,  the  proposed  algorithm  was 
assessed  by  a  series  of  synthetic  images  or  digital  phantoms. 
Typically,  a  thoracic  CT  image  together  with  the  contour  was 
distorted  with  the  intentionally  introduced  deformation,  and 
then  the  contour  was  propagated  onto  the  distorted  image. 


A  quantitative  comparison  was  carried  out.  The  mean  and 
maximum  separation  between  the  gold  standard  and  the  map¬ 
ped  contours  were  found  to  be  1.0  mm  and  1.5  mm,  respec¬ 
tively.  Figure  5  shows  one  example  of  digital  phantom 
experiments. 

The  performance  of  the  proposed  algorithm  was  further 
evaluated  by  the  bidirectional  mapping  calculation  outlined 
in  Methods  and  Materials.  A  template  contour  at  phase  1 
was  first  mapped  to  phases  3  and  6.  The  mapped  contours 
were  then  treated  as  the  “starting  contours”  and  mapped 
back  to  phase  1 .  The  two  back-mapped  contours  were  com¬ 
pared  with  the  original  template  contour.  The  displacement 
of  each  back-mapped  vertex  point  relative  to  its  original  loca¬ 
tions  was  computed,  and  a  mean  value  of  0.8  mm  was  found 
for  the  bidirectional  mapping  between  phases  1  and  3  and  1.8 
mm  between  phases  1  and  6.  The  larger  displacement  in  the 
latter  situation  was  due  to  the  fact  that,  computationally,  it  is 
more  difficult  to  map  between  two  opposite  phases,  such  as 
inhale  and  exhale  phases,  owing  to  larger  organ  deforma¬ 
tions.  Overall,  the  observed  displacement  is  comparable  to 
the  pixel  size,  indicating  that  the  mapping  is  accurate  and 
robust. 


Thoracic  patient  study  results 

Figure  6  shows  the  contour  mapping  results  for  the  first 
clinical  case.  The  results  are  presented  in  axial,  coronal, 
and  sagittal  planes  for  phases  2  (Fig.  6a-c),  6  (Fig.  6d— f),  8 
(Fig.  6g — i),  and  10  (Fig.  6j-l).  For  phases  2  and  10,  which 
are  immediately  adjacent  to  the  inhale  phase,  the  deformation 
is  relatively  small  and  the  mapped  contours  conform  to  the 
ROI  boundary  very  well.  This  represents  the  “easy”  map¬ 
ping  situation  and  is  consistent  with  the  analysis  presented 
above.  The  average  error  was  less  than  1.5  mm.  For  a  “re¬ 
mote”  phase,  such  as  phase  6  shown  in  Fig.  6d-f,  more 


Fig.  5.  Synthetic  image  and  overlaid  contours.  The  original  contour 
is  depicted  in  green,  gold  standard  contour  in  blue,  and  the  mapped 
contour  in  red. 


604 


I.  J.  Radiation  Oncology  •  Biology  •  Physics 


Volume  70,  Number  2,  2008 


Fig.  6.  Computed  tomographic  images  and  mapped  contours  for  thoracic  patient  1.  Displayed  are  selected  phases.  From 
the  top  row  to  bottom,  phases  2,  6,  8,  and  10  are  presented,  respectively.  For  each  phase,  transverse,  coronal,  and  sagittal 
views  are  shown  from  left  to  right. 


iterations  were  entailed  to  find  the  optimal  solution,  and  the 
resultant  contours  tend  to  be  worse  as  compared  with  those 
phases  adjacent  to  phase  1.  According  to  the  bidirectional 
mapping,  the  average  mapping  error  for  phase  6  was  esti¬ 
mated  to  be  less  than  3  mm.  The  mapped  GTV  contours  (in 
red)  together  with  manual  contours  (in  blue)  by  a  physician 


for  phases  1,4,  8,  and  10  in  the  study  of  patient  4  are  shown 
in  Fig.  7  (parts  a,  b,  d,  and  e,  respectively).  The  template 
phase  (phase  6)  with  the  template  manual  contour 
(in  green)  is  shown  in  Fig.  7c.  In  addition,  the  template  man¬ 
ual  contour  from  this  phase  was  overlaid  on  all  the  displayed 
phases.  For  phases  4  and  8  the  defonnation  was  relatively 
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Fig.  7.  Axial  view  of  computed  tomographic  images  with  gross  tumor  volume  contours  for  the  fourth  thoracic  patient,  (a), 
(b),  (c),  (d),  and  (e)  correspond  to  phases  1,4,6  (template  phase),  8,  and  10,  respectively.  The  green  curves  are  the  manually 
outlined  template  contour  from  phase  6,  and  the  red  curves  represent  the  contours  after  warping.  The  manual  contours  (in 
blue)  by  a  physician  on  individual  phases  were  also  displayed. 


small  (the  manual  contour  was  delineated  on  phase  6),  and 
fewer  iterations  were  needed  to  find  the  optimal  bands  on 
the  target  images.  For  phases  1  and  10,  whereby  deformation 
was  significant  in  the  ROIs  although  more  computing  load 
was  necessary,  a  good  result  was  still  achieved  with  our  nar¬ 
row-band  technique.  Comparisons  between  the  mapped  con¬ 
tours  and  the  manually  segmented  contours  by  physicians  for 
these  patients  were  also  performed,  and  results  revealed 
a  similar  level  of  accuracy  (maximum  and  mean  values  of 
the  discrepancy  between  the  two  sets  of  contours  are  2.8 
mm  and  -0.9  mm,  respectively). 

As  a  useful  application  of  the  proposed  technique,  in  Fig.  8 
we  present  the  mean  and  maximum  lung  displacements  of 
contour  vortices  for  each  breathing  phase  relative  to  their  lo¬ 
cations  on  the  template  phase.  As  seen  in  Fig.  8,  the  overall 
behavior  of  the  mean  and  maximum  displacements  is  consis¬ 
tent  with  our  intuitive  expectation.  For  cases  1  and  2,  the  in¬ 
hale  phase  (phase  1)  was  manually  segmented,  thus  the 
displacement  for  that  phase  is  zero.  For  other  phases,  both 
mean  and  maximum  displacement  values  vary  with  the 
breathing  phase  and  reach  their  maxima  at  the  opposite 
phase.  For  case  3  the  exhale  phase  was  manually  segmented, 
and  the  behavior  was  thus  opposite  to  cases  1  and  2.  In  gen¬ 
eral,  an  average  displacement  of  approximately  3  mm  was 
found  for  inhale  and  exhale  phases.  A  slight  digression  is  no¬ 
ticed  in  phase  7  of  patient  1 ,  which  may  be  caused  by  4D-CT 
binning  artifacts.  This  type  of  data  is  particularly  useful  in 
determining  the  patient-specific  tumor  margin  to  account 
for  breathing  motion  of  the  tumor  target. 


Contour  propagation  in  a  head-and-neck  case 

The  results  of  contour  mapping  for  the  head-and-neck  case 
are  summarized  in  Fig.  9.  Figure  9a  shows  the  planning  CT 
along  with  manually  delineated  contours,  and  Fig.  9b  dis¬ 
plays  the  mapped  contours  of  the  body,  mandible,  and 
GTV  on  CBCT.  For  body  and  mandible  a  simple  rigid  map¬ 
ping  is  enough  to  achieve  high  accuracy.  For  the  GTV,  how¬ 
ever,  the  proposed  deformable  registration  model  was 
necessary  to  adequately  propagate  the  contour.  A  visual 
inspection  of  the  propagated  contours  suggests  that  the  map¬ 
ping  is  clinically  acceptable. 

DISCUSSION 

Four-dimensional  CT  image  segmentation  represents 
a  necessary  step  in  constructing  a  4D  patient  model  and  com¬ 
puting  the  accumulated  dose  in  4D  radiotherapy.  A  natural 
way  to  tackle  the  problem  is  to  auto-map  the  manually  delin¬ 
eated  contours  on  one  of  the  phases  to  the  remaining  phases. 
In  this  work,  a  regional  computing  algorithm  was  introduced 
to  deal  with  the  issue.  The  approach  relies  on  the  assumption 
that  a  narrow  band  surrounding  the  manually  segmented  con¬ 
tour  can  capture  sufficient  information  to  drive  the  finding  of 
its  counterparts  in  other  phases  of  the  4D-CT.  Obviously,  this 
assumption  is  valid  when  the  band  is  sufficiently  wide  so  that 
a  large  number  of  voxels  are  involved  in  the  registration  cal¬ 
culation.  As  demonstrated  by  the  presented  data,  the  registra¬ 
tion  and  the  mapping  are  reliable  when  the  bandwidth  is 
larger  than  4  mm.  Computationally,  the  proposed  approach 
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Fig.  8.  Displacement  of  region  of  interest  boundary  points  as  a  function  of  respiration  phase  for  three  thoracic  patients,  (a) 
Mean  displacement  vs.  phase,  (b)  Maximum  displacement  vs.  phase.  4D  CT  =  four-dimensional  computed  tomography. 


resides  between  a  deformable  model-based  mapping  and 
a  surface  model-based  ROI  contour  mapping. 

The  success  of  the  image  content-based  approaches,  such 
as  the  proposed  narrow-band  approach  or  conventional  de¬ 
formable  image  registration,  arises  from  the  fact  that  they 
fully  utilize  the  inherent  image  features  of  the  two  input  im¬ 
ages.  The  narrow  band-based  technique  is  particularly  attrac¬ 
tive  because  it  takes  advantages  of  the  useful  features  of  both 
image  content-based  technique  and  the  regional  surface- 
based  model.  In  a  sense,  it  is  a  hybrid  approach  of  the  two  dis¬ 
tinct  types  of  algorithms.  The  narrow-band  approach  utilizes 
the  imaging  features  surrounding  the  ROI  to  guide  the  search 
of  the  optimal  mapped  contours  while  considering  the  shape 
integrity  of  the  ROI  surface.  It  eliminates  the  need  for  a  global 


registration  of  the  input  images  and  thus  greatly  increases  the 
computational  efficiency. 

Application  of  the  proposed  contour  mapping  technique  to 
five  clinical  cases  indicates  that  the  technique  is  accurate  and 
computationally  efficient.  A  common  problem  in  image 
segmentation  and  contour  mapping  studies  is  the  lack  of 
quantitative  validation.  In  the  studies  of  Lu  et  al.  (13)  and 
Schriebmann  et  al.  (14),  for  example,  the  accuracy  of 
a  deformable  model-based  contour  mapping  technique  was 
evaluated  purely  on  the  basis  of  visual  inspection.  Although 
it  is  a  convenient  way  for  rapid  assessment  of  a  segmentation 
calculation,  especially  in  a  case  in  which  the  “ground  truth” 
contours  do  not  exist,  the  method  falls  short  in  quantization. 
The  same  approach  was  used  in  many  other  previous 


Fig.  9.  Contour  propagation  in  a  head-and-neck  case,  (a)  Planning  computed  tomography  with  manually  outlined  template 
contours  (in  blue)  for  body,  mandible,  and  gross  tumor  volume,  (b)  Cone-beam  computed  tomography  along  with  contours 
after  warping  (in  red)  for  the  corresponding  structures. 
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investigations  (1,5, 14,  30).  In  this  study,  a  bidirectional  con¬ 
tour  mapping  was  proposed  to  examine  the  reliability  and  ro¬ 
bustness  of  a  contour  mapping  technique.  This  method 
provides  a  useful  test  in  assessing  the  success  of  a  contour 
propagation  algorithm.  We  would  like  to  point  out  that  the  bi¬ 
directional  mapping  technique  introduced  in  this  work  is 
a  necessary  (but  not  sufficient)  test.  In  a  rare  but  possible  sit¬ 
uation,  the  bidirectional  mapping  may  not  be  able  to  find  that 
an  error  occurred  in  the  narrow-band  mapping  process.  A  vi¬ 
sual  inspection  of  the  mapped  result  may  help  in  this  situa¬ 
tion.  On  the  basis  of  the  bidirectional  mapping  experiments 
and  visual  inspection  for  the  patient  studies,  we  conclude 
that  the  proposed  approach  can  perform  very  well  even  in 
the  presence  of  significant  deformations. 

In  our  calculation,  we  observed  that  the  regular  grid  of 
BSpline  control  points  could  be  mapped  to  a  region  outside 
the  narrow  band.  Although  it  seems  that  this  does  not  directly 
affect  the  accuracy  of  the  method,  it  may  prolong  the  calcu¬ 
lation  by  computing  the  displacements  in  regions  where  met¬ 
ric  information  is  irrelevant.  Setups  have  been  proposed  to 
adapt  the  splines  control  mesh  to  regions  where  deformation 
is  found  to  be  significant  (31),  and  the  extension  of  the 
method  would  allow  us  to  use  the  BSpline  control  points  de¬ 
fined  only  in  the  regions  within  the  narrow  band.  Implemen¬ 
tation  of  this  type  of  technique  should  further  reduce  the 
computation  time  required  to  find  the  optimal  solution. 

Although  there  are  numerous  deformable  algorithms,  in¬ 
cluding,  for  example,  the  elastic  model  (32-34),  viscous  fluid 
model  (35),  optical  flow  model  (5,30,36),  finite  element 
model  (33,  37),  and  radial  basis  function  models  such  as 
the  basis  spline  model  (28,  38,  39)  and  thin  plate  spline  model 
(40-43),  a  truly  robust  tool  suitable  for  routine  clinical  appli¬ 
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cations  is  yet  to  be  developed.  Each  of  these  approaches  has 
its  pros  and  cons.  The  deformable  calculation  can  be  greatly 
facilitated  if  some  a  priori  system  information  can  be  incor¬ 
porated.  Along  this  line,  the  homologous  correspondence  of 
the  bony  structure  in  two  input  images  has  been  incorporated 
in  thin  plate  spline  method,  and  remarkable  improvement  has 
resulted  (44).  The  narrow  band-generated  ROI  contour  cor¬ 
respondence  could  also  be  used  as  prior  knowledge  to 
improve  a  deformable  registration.  This  work  is  still  in  prog¬ 
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CONCLUSIONS 

In  this  work  we  have  developed  a  regional  deformable 
registration-based  method  to  auto-propagate  contours  for 
4D  radiotherapy.  The  central  idea  is  that  a  narrow  band 
encompassing  an  ROI  surface  carries  the  neighborhood  in¬ 
formation  of  the  ROI  surface  and  can  be  used  to  establish 
a  reliable  association  between  the  ROIs  in  two  phase-specific 
image  sets.  Different  from  other  type  of  regional  algorithms, 
such  as  surface  mapping,  the  method  uses  the  image  features 
captured  in  a  band  to  guide  the  search  for  the  optimal  contour 
mapping.  Compared  with  conventional  deformable  image 
registration-based  approaches,  a  great  reduction  in  computa¬ 
tional  burden  and  a  large  capture  radius  in  optimization  space 
result.  Our  study  demonstrated  that  the  information  contained 
in  the  boundary  region  can  be  used  to  guide  the  contour  map¬ 
ping  in  all  the  testing  cases  presented  in  this  article.  The  pro¬ 
posed  regional  model  decreases  the  workload  involved  in 
4D-CT  ROI  segmentation  and  provides  a  valuable  tool  for 
the  efficient  use  of  available  spatial-temporal  information 
for  4D  simulation  and  treatment  planning. 
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ABSTRACT 

Purpose:  The  purpose  of  this  work  is  to  develop  a  novel  registration  strategy  to 
automatically  map  the  rectal  contours  from  planning  CT  (pCT)  to  cone  beam  CT  (CBCT) 
for  adaptive  radiotherapy. 

Methods  and  Materials:  The  rectal  contours  were  manually  delineated  in  the  pCT.  A 
narrow  band  with  the  delineated  contours  as  its  interior  surface  was  then  constructed. 
The  correspondence  contours  in  the  CBCT  was  found  by  using  a  feature-based 

registration  algorithm,  which  consists  of  two  steps:  (1)  automatically  searching  for 

control  points  in  the  pCT  and  CBCT  based  on  the  feature  of  the  surrounding  tissue  and 
matching  the  homologous  control  points  using  the  Scale  Invariance  Feature 

Transformation  (SIFT);  (2)  using  the  control  points  for  a  Thin  Plate  Spline  (TPS) 

transformation  to  warp  the  narrow  band  and  finding  the  corresponding  contours  from 
pCT  to  CBCT. 

Results:  A  robust  rectal  contour  propagation  method  has  been  developed.  It  was  able  to 
correctly  identify  sufficient  number  of  homologous  control  points  in  the  two  input  images 
and  enabled  accurate  warping  of  the  narrow  band.  Digital  phantom  study  indicated  that 
an  average  accuracy  of  1.2  mm  is  readily  achievable.  For  clinical  cases,  the  method  also 
yielded  satisfactory  results  even  when  there  were  significant  rectal  content  changes 
between  the  pCT  and  CBCT  scans. 

Conclusion:  Exclusion  of  the  volume  inside  the  rectum  and  efficient  detection  of  image 
features  are  two  key  factors  for  successful  rectal  contour  mapping.  The  proposed 
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technique  provides  a  powerful  tool  for  adaptive  radiotherapy  of  prostate,  rectal  and 
gynecological  cancers  in  the  future. 


Key  words:  Image-guided  radiotherapy  (IGRT),  Adaptive  radiotherapy,  Deformable 
registration,  Contour  mapping,  Scale  Invariance  Feature  Transformation  (SIFT). 
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INTRODUCTION 

Patients  treated  with  radiotherapy  for  cancers  such  as  prostate,  rectal,  and 
gynecological  cancers  experience  significant  day-to-day  changes  in  their  rectal  volumes 
due  to  motion,  distention,  and  filling.  Due  to  variations  in  the  image  content,  an  exact 
correspondence  between  two  image  sets  acquired  at  different  time  points  may  not  exist. 
Thus  any  deformable  model  relying  on  the  use  of  infonnation  contained  in  the  entire 
image  may  not  be  adequate  in  dealing  with  these  patients.  The  artifacts-induced  disjoint 
between  the  images  also  makes  the  auto-propagation  of  contours  outlined  in  one  set  of 
images  to  another  highly  difficult  with  conventional  strategies.  With  continued 
enthusiasm  for  adaptive  radiotherapy,  the  ability  to  reliably  and  efficiently  map  the 
rectum  outlined  in  the  pCT  to  the  on-treatment  CBCT  images  now  becomes  a  bottleneck 
and  needs  to  be  resolved  in  order  for  many  patients  with  cancer  within  the  pelvis  to 
benefit  from  the  novel  adaptive  re-planning  strategy  (1,2). 

The  issue  of  rectal  motion  and  deformation  in  conformal  radiation  therapy  is 
described  in  various  publications.  Lee  et  al  (3)  evaluated  the  CBCT  as  a  tool  to  quantify 
the  accuracy  and  precision  of  a  simulated  IMRT  treatment  delivery  model  for  rectal 
cancer  when  rectal  motion  due  to  filling  and  deformation  was  taken  into  account.  The 
mean  deformation  variation  of  0.71  and  0.94  cm  in  the  LAT  and  AP  directions  was 
reported.  Foskey  et  al  (4)  shrank  the  rectal  gas  region  to  a  virtual  point  in  order  to  make 
the  correspondence  of  the  rectal  volumes  in  two  sets  of  images.  Similar  to  that  reported 
by  Schreibmann  et  al  (5),  Gao  et  al  (6)  used  an  automatic  image  intensity  modification 
procedure  to  create  artificial  gas  pockets  in  the  pCT  images.  The  major  drawbacks  of 
these  types  of  approaches  are  the  artificial  introduction  of  image  features  within  the  rectal 
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volume  and  the  potentially  inaccurate  association  of  the  artificial  image  features.  As  a 
consequence,  the  concordance  between  the  rectal  volumes  after  deformable  registration 
and  the  manually  segmented  rectum  was  found  to  be  less  than  80%. 

In  this  work,  we  propose  to  use  the  image  information  in  the  neighborhood 
outside  the  rectal  wall  as  the  driving  force  to  guide  the  rectal  contour  propagation  from 
the  pCT  to  CBCT.  Because  the  content  in  the  region  outside  the  rectal  wall  should  be 
conserved,  regardless  of  any  changes  in  the  rectal  filling  and  distension,  this  strategy 
seems  to  be  physically  sensible.  Coupled  with  a  powerful  feature-based  deformable 
registration  model,  which  identifies  homologous  tissue  features  shared  by  the  pCT  and 
CBCT  images,  the  novel  approach  captures  the  key  issues  of  the  system  and  provides  a 
natural  solution  to  the  above  stated  problem.  Application  of  the  proposed  algorithm  to  a 
number  of  digital  phantoms  and  clinical  cases  demonstrates  that  the  technique  is  accurate 
and  robust  and  may  be  useful  for  future  adaptive  therapy  planning. 


MATERIALS  AND  METHOD 

Software  platform 

The  proposed  contour  mapping  algorithm  was  implemented  using  the  Insight 
Toolkit  (7,  8)  and  the  Visualization  Toolkit  (VTK)  (9),  which  are  open  source  cross¬ 
platform  C++  software  toolkits  sponsored  by  the  National  Library  of  Medicine  (NLM). 
They  are  freely  available  for  research  purposes  (http://www.itk.org  for  ITK  and 
http ://public .kitware .com/YTK/  for  VTK).  ITK  provides  various  basic  algorithms  to 
perform  registration  and  segmentation  for  medical  images.  The  programs  contained  in 
ITK  are  highly  extendable,  making  it  an  ideal  platform  for  development  of  image 
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registration  and  processing  techniques.  VTK  is  primarily  used  for  image  visualization 
(including  contours). 

Narrow  band  representation  of  the  rectal  wall 

Inconsistency  in  rectal  contents  between  two  input  image  sets  could  severely 
reduce  the  perfonnance  of  a  deformable  registration  algorithm.  Co-registering  an  empty 
rectum  without  bowel  gas  to  a  rectum  filled  with  bowel  gas  using  any  deformable  model 
could  be  problematic,  for  example.  A  natural  strategy  is  to  exclude  the  volume  inside  the 
rectal  wall.  In  practice,  the  template  rectal  contour  in  the  pCT  image  has  been  manually 
contoured  as  a  part  of  the  routine  treatment  planning  process,  thus  making  it  a 
straightforward  matter  to  exclude  the  volume  inside  the  rectal  wall.  Figure  1  shows  the 
proposed  contour  mapping  process.  After  manual  segmentation  in  the  pCT,  a  narrow 
band  as  sketched  in  Fig.  2  is  constructed  with  the  manually  segmented  rectum 
representing  the  inner  surface  of  the  band.  On  an  axial  slice,  the  contour  has  a  polygon 
shape  and  the  vertices  of  the  polygon  form  the  basis  for  constructing  the  narrow  band. 
The  distance  between  the  neighboring  vertices  on  the  contour  is  typically  2-10  mm 
depending  on  the  shape  of  the  contour.  In  generating  the  narrow  band,  we  first  create 
cubes  with  side  length  of  d  for  each  vertex,  as  depicted  by  points  A  and  B  in  figure  2b.  In 
order  to  obtain  a  smooth  band,  between  A  and  B  three  more  cubes,  cornered  at  points  C, 
D,  and  E,  are  inserted.  Point  C  is  chosen  to  be  the  middle  point  between  A  and  B.  Point 
D  is  the  point  between  A  and  C,  and  point  E  is  the  point  between  B  and  C.  More 
interpolated  vertex  points  can  be  similarly  introduced  when  needed.  The  blue  area  in  Fig. 


2a  stands  for  the  narrow  band. 
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The  narrow  band  in  our  approach  is  used  as  a  compact  representation  of  the  rectal 
surface.  As  will  be  detailed  in  the  next  subsection,  a  feature-based  deformable 
registration  algorithm  is  employed  to  find  the  correspondence  of  the  band  in  the  CBCT 
images.  Upon  successful  registration,  the  defonnation  field  is  utilized  to  propagate  the 
pCT  contour  to  the  CBCT.  Because  only  the  image  features  outside  the  rectum  is  used,  a 
narrow  band  shown  in  Fig.  2  pennits  us  to  take  advantage  of  the  regional  calculation 
algorithm  yet  avoiding  the  nuisance  of  rectum/bladder  filling. 

Feature-based  warping  of  the  narrow  band 

As  illustrated  in  figure  1,  the  process  of  contour  mapping  is  to  warp  the  narrow 
band  constructed  above  in  such  a  way  that  its  best  match  in  the  CBCT  images  is  found. 
Mathematically,  this  constitutes  an  optimization  problem,  in  which  a  group  of 
transfonnation  parameters  that  transform  the  points  within  the  band  in  pCT  to  their 
corresponding  points  in  CBCT.  The  input  to  the  contour  mapping  software  includes  the 
narrow  band  and  the  CBCT  images,  which  are  described  by  the  image  intensity 
distributions  /a(x)  and  /b(x),  respectively. 

To  find  the  transfonnation  matrix,  T(x),  that  maps  an  arbitrary  point  in  the  band  to 
the  corresponding  point  in  the  CBCT  images  (or  vice  versa),  a  Thin  Plate  Spline  (TPS) 
deformable  model  is  employed.  But  other  models  should  also  be  applicable  to  model  the 
deformation  of  the  band.  Currently,  the  TPS  method  still  needs  manual  placement  of 
control  points  and  this  work  automates  the  control  point  selection  by  using  the  SIFT 
tissue  feature  searching  (see  next  subsection  for  details).  Roughly,  300  control  points  are 
selected  based  on  the  prominent  tissue  features. 
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The  detailed  description  of  the  TPS  transformation  can  be  found  in  Ref.  (10-13). 
Briefly,  a  weighting  vector  W=(w\,  w2. . .  wn)  and  the  coefficients  a\,  aa,  av  are  computed 
from  a  series  of  matrices  which  are  constructed  using  n  pairs  of  SIFT-selected  control 
points  in  the  pCT  image  (x,,  Vi)  and  in  the  CBCT  image  («,,  v,),  respectively.  The 
function  transforming  a  voxel  in  the  pCT  to  a  new  coordinate  in  the  CBCT  is  obtained 
from 

n  ,  v 

f{ii,v')  =  ax  +auu  +  avv  +  YJwiU\\pi-(u,v)  |),  (1) 

i=0 


where  p  is  the  matrix  of  the  control  points  coordinates  in  the  pCT  and  U  is  a  basis 
function  to  measure  the  distance.  Some  major  steps  of  the  TPS  calculation  are: 

( 1 )  Assuming  Px  =  (xt , y1 ) ,  P2  =  (x2,y2) Pn  =  (xn ,  v„ )  are  the  n  control  points 
in  the  pCT  images.  The  distance  between  point  i  and  j  is  given  by  r..  =  \  P:  -  R  | . 
Define  matrices 


1  x1 
1  x2 


y  i 


1  Xn  yn 


(2) 


0  U  (r12)  •••  U(rJ 

U(r2 1)  0  -  U(r2n) 

U(rnl)  U(rn2 )  -  0 


(3) 


and 


K  P 


(4) 


where  O  is  a  3  x3  matrix  of  zeros  and  U  is  a  basic  function  U (r)  =  r  log  r2 . 
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(2)  Letting  Ql  =(w1,v1),  Q2  =(u2,v2),  ...,  Qn  ={un,vn)  be  n  corresponding 
control  points  in  the  CBCT  image.  We  construct  matrices 


V  = 


v,  v2 


(5) 


Y  =  (V  |  0  0  Of. 


(6) 


The  weighting  vector  W  =  (w1,w2,  -,wn)  and  the  coefficients  a i,  au>  and  av  can  be 
computed  by  the  equation 

L~'Y  =  (W  |  a,  au  avf .  (7) 

(3)  Using  the  elements  of  L'lY  to  define  a  function  fin v')  everywhere  as  given  in 
Eq.  (1).  This  function  transfonns  a  voxel  in  the  pCT  volume  to  a  new  coordinate  in  the 
CBCT  image.  Upon  successful  registration,  the  deformation  field  is  extracted  and 
utilized  to  transform  the  manual  rectal  contours  to  CBCT. 


Scale-Invariant  Features  Transformation  (SIFT) 

The  feature-based  deformable  registration  is  an  essential  part  of  the  proposed 
contour  mapping  process.  While  the  TPS  deformable  registration  is  relative  simple  and 
doesn’t  require  iterations  and  intensive  calculation  for  each  individual  voxel,  it  relies  on 
the  use  of  homologous  control  points  in  the  two  input  image  sets  to  be  co-registered.  In 
reality,  the  interactive  identification  of  the  control  points  in  both  images  is  tedious, 
difficult,  and  often  a  source  of  error.  Here  we  automate  the  control  point  selection  by 
using  the  SIFT-based  tissue  feature  searching.  Because  of  the  efficient  use  of  a  priori 
system  knowledge,  the  approach  greatly  enhances  the  robustness  of  the  narrow  band 
warping  algorithm. 
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The  SIFT  method  was  introduced  by  Lowe  (14)  to  characterize  the  local  tissue 
features.  The  method  utilizes  both  image  intensity  and  local  gradient  information  to 
characterize  the  neighborhood  property  of  a  point.  The  algorithm  includes  scale-space 
extrema  detection,  control  point  localization,  orientation  assignment  and  control  point 
descriptor.  In  2D  case,  for  example,  the  method  uses  the  orientation  histograms  of  the 
four  quadrants  surrounding  a  point  (containing  64  pixels)  to  characterize  the  inherent 
tissue  feature  of  the  point  (see  Fig.  3).  To  obtain  the  histogram  for  a  quadrant,  as 
illustrated  in  Fig.  3,  the  gradient  of  each  of  the  16  pixels  in  a  quadrant  is  computed.  An 
eight-bin  histogram,  with  first  bin  representing  the  number  of  pixels  whose  gradients  fall 
between  0°  and  45  °,  and  so  forth,  is  then  constructed.  For  illustration,  the  histogram  of 
each  of  the  four  quadrants  is  displayed  schematically  in  the  right  panel  of  Fig.  3  as  an 
eight-vector  plot.  In  total,  32  vectors  are  calculated  in  2D  case.  In  extending  the  SIFT 
method  from  2D  to  3D,  total  of  192  vectors  are  needed.  These  vectors  represent  the  local 
feature  and  serve  as  a  signature  of  the  point.  The  SIFT  descriptor  is  considered  as  one  of 
the  most  effective  descriptors  currently  available(15,  16). 

Theoretically,  the  SIFT  descriptor  can  be  computed  for  each  voxel  in  an  image. 
However,  this  is  computationally  expensive.  The  commonly  used  sampling  strategy  is 
to  compute  the  descriptor  every  2  to  3  voxels  in  x,  y  and  z  directions.  After  the  SIFT 
descriptors  are  computed  in  both  input  images,  the  points  having  the  most  similar  SIFT 
descriptors  in  the  two  images  are  then  identified.  For  a  given  point,  indexed  by  n,  in  the 
pCT  image,  the  least-squares  difference  of  the  SIFT  descriptor  of  the  point  and  that  of  a 
potential  association  point  n  in  the  CBCT,  SnM-,  is  first  computed  according  to 

V  «=i 
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where  /  represents  the  pixel  intensity,  a  index  the  bins  of  the  SIFT  histogram  of  a  point 
and  the  summation  over  a  runs  from  1  to  32  for  2D  case,  and  192  for  3D  case.  After  Sn>n' 
is  computed  for  all  n  in  the  CBCT,  two  points  n[  and  n'2  that  have  the  least  histogram 

difference  with  point  n  are  identified.  If  the  ratio  of  these  two  values  is  less  than  80%, 
the  point  that  has  the  lowest  S  value  is  chosen  tentatively  as  the  correspondence  of  the 
point  n,  otherwise,  no  association  is  made  for  the  point.  Note  that  the  criterion  of  80% 
here  is  an  empirical  value.  If  the  value  is  too  large,  the  number  of  false  association 
increases.  Conversely,  many  true  associations  may  be  missed.  In  general,  this  criterion 
should  be  determined  on  an  organ  specific  basis.  For  lung,  for  example,  it  was  found  that 
a  criterion  of  50%  is  adequate  to  find  sufficient  number  of  associations. 

To  further  increase  the  accuracy  of  control  point  association,  a  bi-directional 
mapping  strategy  is  developed  based  on  the  fact  that  if  a  point  in  the  pCT  is  mapped 
correctly  to  the  CBCT,  it  will  be  default  to  be  mapped  back  to  the  original  control  point 
in  the  pCT  when  an  inverse  map  is  applied  to  the  corresponding  control  point  in  the 
CBCT.  Therefore,  after  the  original  association  of  control  points  as  described  above,  the 
mapped  points  in  CBCT  is  inversely  co-registered  to  the  pCT.  If  the  correspondence  still 
exists,  the  control  point  pair  is  labeled  a  match.  Otherwise,  they  are  considered  as  a 
mismatch  and  deleted  from  the  list  of  correspondence  points.  Upon  the  association  of  the 
points,  the  associated  points  are  employed  as  control  points  for  TPS  as  described  above. 

Evaluation  of  the  models  using  digital  phantom  and  existing  patient  data 

The  performance  of  the  above  model  is  evaluated  by  a  number  of  2D  digital 
phantoms  and  archived  clinical  cases.  In  the  digital  phantom  experiments,  two 
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deformations  are  introduced.  A  virtue  of  this  approach  is  that  the  “ground  truth” 
solutions  exist  and  the  transformation  matrices  are  known,  thus  making  the  evaluation 
straightforward.  The  mathematical  transfonnations  used  to  deform  the  phantom  is 
generated  using  a  formula  (17): 

x'(x,y')  =  (I  +  bco'imq)x  (9) 

y'(x,y)=  O  +  bcosmq^y  (10) 

Here  q  =  tan  1  - .  Two  parameters,  m  and  b,  are  used  to  characterize  a  deformation. 

X 

Generally,  they  describe  the  complexity  and  magnitude  of  a  deformation,  respectively. 
The  contour  outlined  in  the  original  image  is  then  mapped  to  the  deformed  image.  The 
accuracy  of  the  contour  mapping  calculation  is  assessed  by  comparing  directly  with  the 
contour  from  the  known  transfonnation  matrix. 

Contour  propagation  from  pCT  to  CBCT  is  studied  by  using  a  prostate  cancer  and 
a  rectal  cancer  case.  The  pCT  is  acquired  with  a  GE  Discovery-ST  CT  scanner  (GE 
Medical  System,  Milwaukee,  WI)  approximately  two  weeks  prior  to  the  initiation  of  the 
radiotherapy.  The  on-treatment  CBCT  images  are  acquired  using  the  Varian  Trilogy™ 
(Varian  Medical  Systems,  Palo  Alto,  CA).  Each  slice  of  pCT  or  CBCT  is  discretized  into 
512  x  512  pixels.  The  images  are  transferred  through  DICOM  to  a  high-performance 
personal  computer  (PC)  with  a  Xeon  (3.6  GHz)  processor  for  image  processing.  The 
manually  outlined  contours  in  the  pCT  images  are  mapped  to  CBCT  images  using  the 
proposed  technique.  In  general,  quantitative  validation  of  a  deformable  registration 
algorithm  for  a  clinical  case  is  difficult  due  to  the  lack  of  the  ground  truth  for  clinical 
testing  cases.  For  the  cases  studied  here,  visual  inspection  method  is  employed  to  assess 
the  success  of  the  proposed  algorithm. 
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RESULTS 


2D  digital  phantom  experiment 

The  proposed  algorithm  is  first  tested  using  a  2D  digital  phantom  (Fig.  4a)  with 
two  intentionally  introduced  deformations  of  the  image  shown  in  Figs.  4b  and  4c, 
respectively.  The  rectal  contour  in  pCT  is  manually  delineated  and  shown  in  Fig.  4a  in 
red.  The  deformation  shown  in  Figs.  4b  and  4c  are  obtained  by  setting  the  parameters  b 
and  c  in  Eqs.  (9)  and  (10)  to  ( b  =  2,  m  =  2)  and  ( b  =  2,  m  =  3),  respectively.  The  green 
contours  in  Figs.  4b  and  4c  represent  the  auto-mapped  contour.  Overall,  the  mapped 
contours  can  capture  the  main  features  of  the  two  dramatic  defonnations,  and  conform 
snugly  to  the  boundary  of  the  rectum  in  both  cases. 

In  obtaining  the  result  shown  in  Fig.  4b,  a  total  of  200  control  points  were 
identified  by  the  bi-directional  SIFT  calculation  as  described  in  method.  For  clarity,  a 
selection  of  the  SIFT-identified  control  point  associations  are  displayed  in  Fig.  5.  The 
total  number  of  control  points  identified  here  are  far  more  than  that  commonly  used  in 
TPS  calculation  (11),  allowing  an  improved  deformable  warping  of  the  narrow  band. 
The  displacement  field  derived  by  using  the  TPS  method  is  shown  in  Fig.  6a.  For 
comparison,  the  known  displacement  field  from  Eqs.  (9)  and  (10)  is  plotted  in  Fig.  6b. 
The  subtraction  between  the  TPS-derived  displacement  field  and  the  known  field  is 
shown  in  Fig.  6c.  It  is  found  that  the  average  deviation  of  the  SIFT-TPS  displacement 


from  the  known  solution  is  less  than  1 .2mm. 
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Clinical  case  study 

The  contour  propagation  study  from  pCT  to  CBCT  for  the  prostate  case  is 
presented  in  Fig.  7.  The  top  row  shows  the  pCT  image  with  manually  delineated 
contours  (green  curves).  The  auto-mapped  contours  overlaid  in  the  CBCT  are  displayed 
in  the  bottom  row  (red  curves).  As  mentioned  in  the  introduction,  the  propagation  of 
rectum  wall  is  often  complicated  by  the  fact  that  the  physical  one-to-one  correspondence 
may  not  exist  due  to  the  addition  or  subtraction  of  some  contents  within  the  rectum. 
Figure  8  exemplifies  this  and  shows  that  the  rectal  filling  at  the  time  of  CBCT  acquisition 
is  quite  different  from  that  of  pCT.  As  can  be  intuitively  conceived,  this  image  content 
change  could  severely  reduce  the  performance  of  a  conventional  deformable 
registration(18-20).  The  narrow  band  approach  described  in  this  work  circumvents  the 
problem  by  excluding  the  rectal  volume  affected  by  the  rectum/bladder  filling.  As  a 
result  in  Fig.  7,  the  mapped  contours  closely  conform  to  the  rectal  wall  change  and  the 
final  contours  are  clearly  clinically  sensible. 

In  practice,  rectal  volume  motion  and  defonnation  can  cause  significant 
uncertainties  pertaining  to  the  adequacy  of  actual  dose  delivered  to  the  gross  tumor 
volume  (GTV)  as  well  as  to  the  surrounding  normal  structures.  This  issue  has  been  a 
major  obstacle  in  the  implementation  of  IMRT  in  rectal  cancer.  In  Fig.  8  the  three  axial 
pCT  and  CBCT  images  of  a  rectal  cancer  patient  acquired  in  an  interval  of  two  weeks  are 
shown.  Significant  target  volume  motion  and  deformation  are  observed  from  Fig.  8.  The 
rectal  volume  in  the  pCT  is  found  to  be  more  than  3  times  that  of  the  rectal  volume  in  the 
CBCT  and  thus  represents  a  challenging  situation  for  any  deformable  model.  The  rectal 
contours  are  manually  drawn  in  the  pCT  and  mapped  to  the  subsequent  CBCT  using  the 
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proposed  method.  The  upper  row  of  Fig.  8  shows  three  axial  slices  of  the  pCT  with 
manually  delineated  contours  (green  curves).  The  results  of  contour  propagation  from 
the  pCT  to  the  CBCT  are  shown  in  the  lower  row  of  Fig.  8  (red  curves). 

DISCUSSION 

In  this  work,  an  effective  feature-based  rectal  contour  mapping  algorithm  has  been 
described.  An  indispensable  step  toward  online  or  offline  adaptive  re -planning  with 
consideration  of  the  patient’s  dose  delivery  history  and  on-treatment  anatomy  is  the 
expedite  organ  segmentation  of  CBCT  images  (7,  21-25).  While  this  task  is,  in  principle, 
achievable  using  deformable  registration  of  the  pCT  and  CBCT  images,  the  accuracy  of 
the  registration  and  therefore  the  contour  mapping,  is  often  adversely  affected  by  the 
presence  of  image  contents  in  one  image  that  do  not  have  correspondence  in  the  other 
image.  The  propagation  of  rectum  wall  is  an  example  of  this.  For  prostate,  rectal,  or 
gynecological  cancer  patients  for  example,  the  presence  and  absence  of  bowel  gas  can 
vary  daily.  Co-registering  an  empty  rectum  without  bowel  gas  to  a  rectum  filled  with 
bowel  gas  (or  vice  versa)  using  any  deformable  model  could  be  problematic  and  large 
errors  could  occur. 

We  describe  a  regional  contour  propagation  algorithm  taking  into  account 
possible  organ  deformation  and  anatomic  changes.  Because  the  narrow  band  contains 
only  the  image  features  outside  the  rectum,  this  method  is  not  affected  by  the  rectum 
filling  changes.  The  proposed  approach  relies  on  the  assumption  that  a  narrow  band 
surrounding  the  manually  segmented  rectal  contour  can  capture  sufficient  information  to 
drive  the  finding  of  its  counterpart  in  the  subsequent  CBCT.  Obviously,  this  assumption 
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is  valid  when  the  band  is  sufficiently  wide  so  that  enough  image  features  are  involved  in 
the  warping  calculation. 

The  use  of  SIFT  descriptor  enhances  our  ability  to  find  the  correspondence  of  the 
narrow  band  because  of  the  effective  utilization  of  image  intensity  and  gradient 
information.  In  contrast  to  the  conventional  intensity-based  image  registration,  which 
only  uses  intensity  infonnation  of  the  voxels,  the  feature-based  registration  extracts 
information  regarding  image  structure,  including  shape,  texture,  etc.  Therefore,  the 
feature-based  image  registration  is  generally  more  effective  in  correctly  identifying 
corresponding  voxels  compared  to  the  intensity-based  image  registration. 

The  proposed  contour  mapping  technique  is  applied  to  digital  phantoms  and 
clinical  cases  and,  in  all  cases,  the  contour  mapping  results  are  found  to  be  clinically 
acceptable.  It  is  important  to  emphasize  that  in  these  test  cases,  the  rectal  deformations 
are  quite  significant  and  thus  present  challenges  to  any  deformable  model  or  contour 
mapping  technique.  It  is  impressive  that  a  simple  approach  with  a  narrow  band  and  SIFT 
descriptor  can  capture  the  main  feature  of  the  rectal  contour  and  help  to  find  the 
correspondence  contours  in  the  CBCT  images. 

In  this  study,  a  bi-directional  SIFT  descriptor  are  employed  to  examine  the 
reliability  and  robustness  of  the  calculations.  The  bi-directional  mapping  further 
enhances  the  degree  of  success  of  a  contour  propagation  algorithm.  It  is  useful  to  note 
that  the  bi-directional  mapping  is  a  necessary  (but  not  sufficient)  test.  In  a  rare  but 
possible  situation,  the  bi-directional  mapping  may  not  be  able  to  find  an  error  occurred  in 
the  contour  mapping  process.  A  visual  inspection  of  the  mapped  result  is  always  helpful 
in  practice. 
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CONCLUSION 

Significant  inter-fractional  patient  setup  uncertainty  and  anatomy  changes  have 
been  reported  in  numerous  studies,  and  are  widely  recognized  as  one  of  the  major 
limiting  factors  for  maximum  exploitation  of  modern  radiation  therapy  techniques  such 
as  IMRT  and  IGRT.  The  advent  of  onboard  volumetric  imaging  devices  promises  to 
improve  the  situation  by  providing  valuable  3D  (or  even  possibly  4D)  geometric  data  of 
the  patient  in  the  treatment  position  and  allows  for  the  adaptive  modification  of  treatment 
plan  during  a  course  of  treatment. 

In  this  work,  an  effective  feature-based  rectal  contour  mapping  algorithm  has  been 
described.  The  method  yielded  satisfactory  mapping  for  both  digital  phantom  and  clinical 
cases.  It  is  impressive  that  the  algorithm  is  able  to  successfully  map  the  contours  from 
pCT  to  CBCT  even  for  some  very  challenging  cases  in  which  the  deformation  and/or 
image  content  change  are  dramatic.  The  two  salient  features  of  the  described  algorithm 
are:  (1)  the  use  of  inherent  tissue  feature  as  a  priori  knowledge  for  deformable 
registration;  and  (2)  limiting  the  ROI  to  exclude  the  volume  inside  the  rectum  and 
focusing  on  the  adjacent  neighborhood  of  the  rectal  contour.  The  algorithm  should  be 
extendable  for  contour  propagation  of  organs  with  similar  features,  such  as  the  bladder 


and  stomach. 
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FIGURE  CAPTIONS 

Figure  1  Overall  process  of  rectal  contour  propagation. 

Figure  2  A  sketch  of  narrow  band,  (a)  a  narrow  band  image  surrounding  a  manually 
segmented  rectal  contour;  (b)  a  narrow  band  construction  is  illustrated  for  two  vortex 
points  A  and  B. 

Figure  3  A  sketch  of  orientation  histogram  in  SIFT  method. 

SIFT — Scale  Invariance  Feature  Transformation 

Figure  4  Rectal  contour  propagation  from  the  pCT  to  two  dramatically  deformed  images, 
(a)  original  contour,  the  red  curve  represents  the  manually  delineated  contour;  (b)  and  (c) 
its  optimal  mapping  in  the  two  defonned  images  (green  curves).  For  comparison,  the 
original  contour  from  the  pCT  is  also  shown  in  the  deformed  images  (red  curves). 

PCT — planning  CT 

Figure  5  Control  points  in  the  2D  contour  mapping 
2D — two  dimentional 

Figure  6  Displacement  fields,  (a)  TPS-derived  displacement  field  for  the  2D  digital 
phantom  study;  (b)  intentionally  introduced  displacement  field;  and  (c)  subtraction  of 
TPS  derived  and  the  known  displacement  fields. 
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TPS — Thin  Plate  Spline 
2D — two  dimentional 

Figure  7  3D  contour  mapping  for  the  rectum  of  a  man  with  prostate  cancer.  The  top  row 
is  the  three  transactions  in  the  planning  CT  image  with  manually  delineated  contours 
(green  contours),  the  bottom  row  is  corresponding  transactions  in  the  CBCT  image  with 
auto-mapped  contour  (red  contours).  The  left  column  is  the  axial  plane,  the  middle 
column  is  the  coronal  plane,  and  the  right  column  is  the  sagittal  plane. 

3D — three  dimentional 

CT — computational  tomography 

CBCT — cone  beam  computational  tomography 

Figure  8  Rectal  contour  mapping  for  a  rectal  cancer  case.  The  top  row  shows  several 
axial  slices  in  the  pCT  image  with  manually  delineated  contours  (green  contours).  The 
bottom  row  is  the  corresponding  slices  in  the  CBCT  image  with  auto-mapped  contours 
(red  contours). 

PCT — planning  CT 

CBCT — cone  beam  computational  tomography 
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ABSTRACT 


The  purpose  of  this  work  is  to  develop  an  effective  technique  to  automatically  propagate 
contours  from  planning  CT  to  cone  beam  CT  (CBCT)  to  facilitate  CBCT  guided  prostate 
adaptive  radiation  therapy.  Different  from  other  disease  sites,  such  as  the  lungs,  the 
contour  mapping  here  is  complicated  by  two  factors:  (i)  the  physical  one-to-one 
correspondence  may  not  exist  due  to  the  insertion  or  removal  of  some  image  contents 
within  the  region  of  interest  (ROI);  and  (ii)  reduced  contrast  to  noise  ratio  of  the  CBCT 
images  due  to  increased  scatter.  To  overcome  the  issues,  we  investigate  a  strategy  of 
excluding  the  regions  with  variable  contents  by  a  careful  design  of  a  narrow  shell 
signifying  the  contour  of  an  ROI.  For  rectum,  for  example,  a  narrow  shell  with  the 
delineated  contours  as  its  interior  surface  was  constructed  to  avoid  the  adverse  influence 
of  the  day-to-day  content  change  inside  the  rectum  on  the  contour  mapping.  The 
corresponding  contours  in  the  CBCT  were  found  by  warping  the  narrow  shell  though  the 
use  of  BSpline  defonnable  model.  Both  digital  phantom  experiments  and  clinical  case 
testing  were  carried  out  to  validate  the  proposed  ROI  mapping  method.  It  was  found  that 
the  approach  was  able  to  reliably  warp  the  constructed  narrow  band  with  an  accuracy 
better  than  1.3  mm.  For  all  five  clinical  cases  enrolled  in  this  study,  the  method  yielded 
satisfactory  results  even  when  there  were  significant  rectal  content  changes  between  the 
planning  CT  and  CBCT  scans.  The  overlapped  area  of  the  auto-mapped  contours  over 
90%  to  the  manually  drawn  contours  is  readily  achievable.  The  proposed  approach 
pennits  to  take  advantage  of  the  regional  calculation  algorithm  yet  avoiding  the  nuisance 
of  rectum/bladder  filling  and  provide  a  useful  tool  for  adaptive  radiotherapy  of  prostate  in 
the  future. 


1.  Introduction 

There  are  increased  interests  in  online/off-line  adaptive  replanning  with  consideration  of 
updated  patient  anatomy  infonnation  derived  from  cone  beam  CT  (CBCT)(Yan  et  al. 
2000,  Court  et  al.  2005,  Ghilezan  et  al.  2005,  Langen  et  al.  2005,  Letourneau  et  al.  2005, 
Oldham  et  al.  2005,  Xing  et  al.  2006,  Yang  et  al.  2007).  An  indispensable  step  toward  the 
realization  of  this  type  of  adaptive  replanning  is  the  segmentation  of  CBCT  images(Gao  et  al. 
2006,  Smeenk  et  al.  2007).  While  this  task  is,  in  principle,  achievable  using  defonnable 
registration  of  the  planning  CT  and  CBCT  images,  the  accuracy  of  the  registration,  and 
therefore  contour  mapping,  is  often  adversely  affected  by  the  presence  of  image  contents  in 
one  image  that  do  not  have  correspondence  in  the  other  image,  such  as  the  presence/absence 
of  bowel  gas  or  bladder  filling  in  the  pelvic  patients(Huang  et  al.  2002).  A  regional  contour 
mapping  algorithm  is  highly  desirable  for  its  high  computational  efficiency  and 
accuracy(Schreibmann  and  Xing  2005,  Chao  et  al.  2007,  Chao  et  al.  2008).  It  also  has  the 
potential  to  improve  the  robustness  of  defonnable  registration  because  the  mapped  contours 
can  be  utilized  as  a  priori  system  infonnation  to  facilitate  the  registration  of  the  two  input 
images. 

Large  interfraction  organ  motion  exists  in  radiation  therapy  of  pelvic  diseases, 
such  as  prostate  and  cervical  cancers(Balter  et  al.  1995,  Yan  and  Lockman  2001, 
Michalski  2003).  Adaptive  therapy  provides  a  viable  option  to  ensure  an  adequate  dose 
coverage  of  the  tumor  target  while  sparing  the  sensitive  structures(Yan  et  al.  1997,  Court 
et  al.  2005,  Wu  et  al.  2008).  Toward  the  general  goal  of  clinical  implementation  of 
adaptive  therapy  of  prostate  cancer,  in  this  work  we  develop  an  effective  regional  contour 
mapping  technique  to  automatically  propagate  rectum  and  prostate  contours  from 
planning  CT  to  CBCT.  Different  from  other  organs,  such  as  the  lungs  and  liver,  the 
contour  mapping  here  is  complicated  by  two  factors:  (i)  the  physical  one-to-one 
correspondence  may  not  exist  due  to  the  insertion  or  removal  of  some  image  contents 
within  the  rectum(Gao  et  al.  2006);  and  (ii)  reduced  contrast  to  noise  ratio  (CNR)  of  the 
CBCT  images  due  to  increased  scatter  in  CBCT  scanning(Li  et  al.  2007).  A  solution 
customized  to  rectum  mapping  is  proposed  to  overcome  the  above  two  limitations  and 


allow  us  to  take  advantage  of  the  regional  calculation  algorithm  yet  avoiding  the  nuisance 
of  rectum/bladder  filling. 


2.  Methods  and  Materials 

Software  platform 

The  narrow  band  contour  mapping  algorithm  was  implemented  using  the  NLM 
Insight  Toolkit  (ITK)(Ibanez  et  al.  2003)  and  the  Visualization  Toolkit  (VTK)(Schroeder 
et  ah),  which  are  open  source  cross-platform  C++  software  toolkits.  They  are  freely 
available  for  research  purposes  (http://www.itk.org  for  ITK  and 
http://public.kitware.com/VTK/  for  VTK).  ITK  provides  various  kinds  of  basic 
algorithms  for  performing  registration  and  segmentation  for  medical  images.  The  fact  that 
programs  contained  in  ITK  are  highly  extendable  makes  it  an  ideal  platform  for  the 
development  of  image  registration  methods.  VTK  is  primarily  used  for  image 
visualization. 

The  contour  propagation  process 

The  rectum  contour  propagation  from  planning  CT  to  CBCT  was  carried  out  in 
the  following  sequences:  (1)  rigid  alignment  of  the  planning  CT  and  CBCT  images  using 
a  rigid  registration;  (2)  construction  of  a  narrow  band  extended  to  a  region  l~2cm  outside 
the  manually  delineated  contours  on  planning  CT;  (3)  warping  of  the  shell  from  the 
planning  CT  to  CBCT.  Upon  successful  warping  of  the  shell,  the  deformation  field  in 
step  (3)  was  utilized  to  transform  the  manually  segmented  contours  to  the  CBCT  images. 
The  central  idea  here  is  to  exclude  the  volume  inside  the  rectum  because  its  content  may 
change  from  day  to  day.  Prostate  contour  can  be  mapped  similarly,  but  the  shell  spans  the 
regions  both  inside  and  outside  the  manually  segmented  prostate.  Figure  1  depicts  the 
flowchart  of  the  overall  contour  mapping  process  using  the  narrow  shell  technique. 

The  manually  delineated  ROI  contours  on  planning  CT  are  referred  to  as  the 
template  contours.  Figure  2  illustrates  the  planning  CT  and  the  template  contours  for  the 
ROIs  such  as  prostate  and  rectum  for  the  one  of  the  patients  studied  in  this  work.  A 


narrow  shell  was  constructed  using  the  method  described  above  for  each  structure. 
Typical  shapes  of  the  narrow  shells  for  rectum  are  illustrated  in  Fig.  3.  The  shell  width 
was  chosen  to  be  10  mm  in  this  study. 

Upon  the  construction  of  shell  structure,  the  shell  is  warped  onto  CBCT  images 
by  using  a  BSpline  model(Lee  et  al.  1996,  Lee  et  al.  1997,  Badea  et  al.  2008)  with  a 
nonnalized  cross  correlation  metric  function  defined  as 


where  lf{x)  and  lf(x')  are  the  intensities  of  the  i-th  voxels  of  images  A  and  B, 
respectively;  x  =  Tx,  where  T  is  the  transformation  matrix.  In  the  narrow  shell  warping 
calculation,  the  input  images  to  the  registration  software  are  the  narrow  shell,  IA(x)  ,  and 

the  CBCT  image,  IB  (x') .  The  narrow  shell  acts  as  a  shape  representation  model  of  an 
anatomic  structure.  The  task  here  is  to  find  the  deformation,  T(x),  that  maps  an  arbitrary 
point  x  on  the  fixed  image  to  the  corresponding  point  x  on  the  moving  image  (or  vice 
versa)  so  that  the  best  possible  match,  as  measured  by  the  registration  metric,  is  achieved. 
The  calculation  proceeds  in  an  iterative  fashion.  The  limited  memory  Broyden-Fletcher- 
Goldfarb-Shannon  algorithm  (L-BFGS)(Liu  and  Nocedal  1989)  was  used  for  the  iterative 
optimization.  The  L-BFGS  optimizer  is  well  known  for  its  superior  perfonnance  in 
dealing  with  high  dimensional  problems  (Schreibmann  and  Xing  2005,  Schreibmann  et 
al.  2006).  We  briefly  describe  the  optimization  algorithm  in  the  following. 

The  BFGS  method  is  derived  from  the  Newton’s  method  in  optimization  which  is  a 
class  of  hill-climbing  optimization  techniques.  It  tries  to  seek  the  stationary  point  of  a 
function,  where  the  gradient  is  0.  Starting  from  a  positive  definitive  approximation  of  the 
inverse  Hessian  Ho  at  xo,  L-BFGS  derives  the  optimization  variables  by  iteratively 
searching  through  the  solution  space.  At  an  iteration  k,  the  calculation  proceeds  as 
follows: 


1 .  Determine  the  descent  direction  pA  =  -H kVf  (xk)  ; 

2.  Line  search  with  a  step  siz eak  =argmin  f(xk  +  cxpk),  where  a  is  the  step  size 

a>  0 

defined  in  the  L-BFGS  software  package; 

3 .  Update  xk+1  =  xk  +  ak  pk  ; 

4.  Compute  Hk+i  with  the  updated  Hk . 

At  each  iteration  a  backtracking  line  search  is  used  in  L-BFGS  to  detennine  the  step  size 
of  movement  to  reach  the  minimum  of/  along  the  ray  xk  +  apk.  During  optimization,  the 
iterative  calculation  continues  until  the  following  stopping  criterion  is  fulfilled: 

ii  y/(xk)  1I2 

max(l,||x,  ||2) 

or  a  pre-set  maximum  number  of  iterations  is  reached.  In  this  study  we  set  s  =  1 0  6  and 
the  iteration  number  to  200. 

Evaluation  of  algorithm  performance 

Evaluation  of  a  contour  propagation  algorithm  is  often  a  difficult  task  because  of 
the  lack  of  the  ground  truth  for  comparison.  A  commonly  used  approach  is  to  visually 
inspect  the  mapped  contours  by  a  clinician.  The  use  of  synthetic  images  (digital  phantoms) 
is  also  useful.  In  this  test,  the  images  and  existing  contours  are  distorted  with  known 
deformation  fields,  which  serves  as  the  ground  truth  for  quantitative  evaluation  of  the 
model  based  contour  propagation.  In  constructing  the  digital  phantom,  the  deformation 
field  used  to  defonn  a  planning  pelvic  CT  image  was  generated  as  follows.  A  B-Spline 
8x8x8  grids  were  overlaid  on  the  CT  image.  Three  random  generators  following  the 
Gaussian  distribution  were  called  for  to  provide  the  displacements  for  the  components  of 
x,  y,  and  z  at  each  grid  point.  The  mean  value  //  of  the  Gaussian  function  were  all 
chosen  to  be  0,  but  the  variances  of  the  three  components  were  set  to  be 


cr  =  u 2  =5.0  mm 

x  y 

a1.  =  2.0  mm. 


Figure  4  (a)  shows  the  deformation  field  overlaid  on  the  original  CT  image.  The  manual 
contour  on  the  original  CT  image  was  warped  with  the  same  defonnation  field.  This 
contour  served  as  the  ground  truth.  A  comparison  of  the  mapped  contour  with  the  ground 
truth  allowed  us  to  quantitatively  assess  the  accuracy  of  the  proposed  approach. 

Patient  Study 

The  image  data  sets  for  five  prostate  cancer  patients  were  used  to  further  evaluate 
the  proposed  technique.  The  CBCT  images  were  acquired  using  a  Trilogy™  system 
(Varian  Medical  Systems,  Palo  Alto,  CA).  All  patient  CBCT  images  were  acquired  with 
90  inA  tube  current  and  120  kV  voltage.  Under  this  image  acquisition  protocol,  the  dose 
was  around  4cGy  for  a  kV  CBCT  scan  (Wen  et  al.  2007).  The  CBCT  image  sets  for  all 
enrolled  patients  were  reconstructed  with  a  2.5  mm  slice  thickness.  The  planning  CT 
scans  were  acquired  with  a  GE  Discovery-ST  CT  scanner  (GE  Medical  System, 
Milwaukee,  WI).  The  size  of  a  CT  slice  was  512  x  512.  The  planning  CT  images  were 
reconstructed  with  a  1.25  mm  slice  thickness  for  all  patients.  All  CT  images  were 
transferred  through  DICOM  to  a  high-performance  personal  computer  (PC)  with  Xeon 
(3.6  GHz)  processor  for  image  processing. 


3.  Results  and  Discussions 

I.  Digital  phantom  study 

The  narrow  shell  constructed  based  on  the  planning  CT  image  and  the  template 
rectal  contour  (Fig.  4(a))  was  warped  to  the  digitally  deformed  images  using  the 
technique  described  above.  For  convenience,  the  deformation  field  used  to  deform  the 
planning  CT  images  is  displayed  in  Fig.  4(a)  as  a  vector  map.  The  rectal  contour  was  then 
extracted  upon  the  successful  warping  of  the  shell.  The  resultant  rectum  contour  is  shown 
in  Figure  4(b)  together  with  the  ground  truth  contour.  The  difference  between  the  two 
sets  of  contours  was  found  to  be  insignificant  by  visual  inspection.  A  quantitative 


comparison  indicated  that  the  mean  discrepancy  between  the  two  sets  of  contours  was 
less  than  1 .3  mm. 

II.  Clinical  case  study 

The  proposed  approach  was  applied  to  five  enrolled  prostate  cancer  patients. 
Figure  5  shows  the  result  of  contour  propagation  for  the  first  clinical  case  in  axial, 
coronal,  and  sagittal  views  of  CBCT  image,  respectively.  For  comparison  purpose,  the 
result  obtained  by  registering  the  whole  images  using  conventional  BSpline  method  was 
also  superimposed  in  Fig.  5  (indicated  by  the  yellow  curves).  We  also  observed  the 
convergence  behaviors  of  the  two  approaches  by  examining  the  metric  function  as  a 
function  of  the  iteration  step.  Nearly  30  iterations  were  needed  in  the  conventional 
registration,  whereas  only  12  iterations  were  necessary  to  obtained  converged  result  in 
narrow  shell  based  approach.  In  this  particular  case,  both  approaches  yielded  acceptable 
results  because  no  substantial  change  in  rectal  filling  in  planning  CT  and  CBCT.  The 
difference  between  two  sets  of  mapping  is  well  within  2  mm.  However,  the  narrow  shell 
based  calculation  was  found  to  be  an  order  of  magnitude  more  efficient  as  measured  by 
the  computational  speed,  in  addition  to  the  greatly  reduced  usage  of  computer  memory. 
Results  for  the  second  patient  are  shown  in  Figure  6.  For  this  patient,  however,  the  rectal 
fillings  in  planning  CT  and  CBCT  were  different  and  the  conventional  BSpline  approach 
lead  to  a  larger  error  (mean  error  ~  8.0  mm)  or  failed  to  accurately  map  the  rectum  from 
planning  CT  to  CBCT.  The  physician’s  manual  contour  was  also  superimposed  on  the 
images  for  visual  comparison  purpose. 

The  mapped  rectum  and  prostate  contours  for  another  case  are  displayed  in  Figure 
7.  Overall,  the  algorithm  performed  well  for  both  rectum  and  prostate.  The  mean  error  of 
the  mapped  contours  was  estimated  to  be  around  2  mm  for  the  rectum,  and  2.5mm  for  the 
prostate.  The  slightly  larger  error  in  prostate  contour  mapping  is  attributed  to  that  fact 
that  the  image  contrast  in  the  prostate  boundary  region  was  not  high.  In  this  case,  a 
careful  examination  and  even  a  manual  modification  of  the  mapped  contour  by  a 
physician  may  be  needed  in  clinical  operation. 

Up  to  this  point,  little  has  been  done  to  register  images  with  the  general  problem 
of  image  content  change  as  studied  in  this  work.  Gao  et  al  and  Fosekey  et  al  (Foskey  et 


al.  2005,  Gao  et  al.  2006)  studied  the  topics  by  using  different  models.  The  essence  of 
their  work  was  to  introduce  “artificial  gas”  to  mimic  the  rectal  image  content  or  shrink 
the  “gas”  region  by  deflation.  While  the  introduction  of  the  artificial  change  could 
improve  the  registration  accuracy  in  some  special  cases,  the  approaches  may  fail  when 
the  change  in  rectum  content  involves  of  solid  substances,  such  as  the  bowel.  A  narrow 
shell  model  seems  to  be  natural  in  addressing  the  issue  of  image  content  inconsistence. 
The  technique  is  conceptually  simple,  easy  to  implement  and  valid  for  a  wide  variety  of 
clinical  situations. 

One  of  the  practical  concerns  is  that  the  relatively  low  quality  of  CBCT  images 
may  influence  the  accuracy  of  image  registration  and  thus  the  contour  mapping.  Paquin  et 
al.  (Paquin  et  al.  2008)  quantitatively  studied  the  influence  of  different  types  of  noises  on 
deformable  registration  and  found  that  the  accuracy  of  image  registration  does  not  depend 
on  the  global  noise  unless  the  noise  reaches  a  certain  threshold  value.  Murphy  et  al. 
(Murphy  et  al.  2008)  also  demonstrated  that  noise  levels  in  cone-beam  CTs  that  might 
reduce  manual  contouring  accuracy  do  not  reduce  image  registration  and  automatic 
contouring  accuracy. 

In  principal,  the  proposed  method  is  also  applicable  to  facilitate  the  contour 
warping  from  planning  CT  to  MV  CBCT  or  images  acquired  using  the  Tomotherapy 
devices.  However,  it  is  important  to  emphasize  that  a  prerequisite  of  any  defonnable 
model  based  approach  is  that  there  are  sufficient  image  features  in  the  neighborhood  of 
the  contours  to  guide  the  warping  of  the  images.  Because  of  the  inherent  low  contrast  of 
the  MV  images,  the  deformable  registration  of  CT  and  MV  CBCT  may  be  challenging. 
Lu  et  al.  have  reported  successful  deformable  registration  of  CT  and  MV  CT  images  (Lu 
et  al.,  Phys.  Med.  Biol.  49(2004)  3067-3087).  It  will  be  an  interesting  subject  of  future 
research  to  examine  whether  a  similar  or  better  accuracy  of  registration  can  be  achieved 
when  a  local  deformable  registration  technique  is  used  for  the  system. 


5.  Summary 


An  automatic  contour  mapping  technique  customized  for  adaptive  radiation  therapy  of 
prostate  cancer  has  been  described.  Deformable  image  registration  has  been  extensively 
studied.  The  narrow  shell  constructed  surrounding  a  ROI  contour  provides  a  compact 
representation  of  the  ROI  and  greatly  facilitates  the  contour  mapping  calculation.  The 
chief  advantage  of  the  proposed  technique  is  that  it  is  much  less  prone  to  the  variation  in 
the  image  contents,  which  occurs  frequently  in  the  pelvic  region  due  to  involuntary 
physiological  process.  Additionally,  the  approach  is  computationally  more  efficient  with 
lower  computer  memory  consumption  and  fast  computational  speed.  Both  phantom  and 
clinical  case  studies  yielded  satisfactory  results  and  suggested  that  the  method  may  prove 
to  be  useful  for  adaptive  radiotherapy  of  prostate  cancer  in  the  future. 
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FIGURE  CAPTIONS 


Figure  1 .  Flowchart  of  the  narrow  shell  contour  mapping  process. 

Figure  2.  Axial  view  (a),  coronal  view  (b),  sagittal  view  (c),  and  volumetric  rendering  (d) 
of  planning  CT  of  the  first  patient.  Template  contours  of  prostate  (magenta)  and  rectum 
(green)  were  superimposed  on  CT  images. 

Figure  3.  Narrow  shell  representation  of  regions  of  interest  (ROIs).  The  curves  represent 
contours  and  the  shaded  regions  sketch  the  narrow  shell  construct.  Illustrated  are  axial 
view  (a),  coronal  view  (b)  and  sagittal  view  (d)  of  the  planning  CT  image  with  a  single¬ 
sided  narrow  shell  surrounding  the  rectum  template  contour. 

Figure  4.  Digital  phantom  experiment,  (a)  original  CT  image  and  manual  contour;  The 
vectors  overlaid  on  the  image  are  the  displacement  field  introduced  to  deform  the  image, 
(b)  digitally  defonned  CT  image  with  ground  truth  contours  and  mapped  contours.  The 
template  contour  after  rigid  registration  of  the  two  images  is  also  displayed. 

Figure  5.  Rectal  contours  on  axial,  coronal,  and  sagittal  planes  obtained  with  the 
proposed  narrow  shell-based  (in  red)  and  conventional  BSpline  approach  (in  yellow)  for 
the  first  clinical  case.  The  contours  after  rigid  alignment  (in  blue)  were  overlaid  on  the 
same  CT  images. 

Figure  6.  Propagated  rectal  contour  (in  red)  overlaid  on  CBCT  images  in  the  second 
patient  case.  The  results  based  on  conventional  approach  (in  yellow)  and  physician’s 
manual  contour  were  also  superimposed  on  CBCT  images. 

Figure  7.  Contour  propagation  results  for  the  fourth  clinical  case,  (a),  (b),  and  (c)  show 
the  rectum  contours  on  axial,  coronal  and  sagittal  views  of  CBCT  image;  (d),  (e),  and  (f) 
display  the  mapped  prostate  contours  on  CBCT  images.  The  manual  contours  were  also 
overlaid  on  CT  images.  The  rectum  contours  are  in  green  and  prostate  in  magenta. 
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